Method and apparatus for simulating physical fields

ABSTRACT

In order to design on-chip interconnect structures in a flexible way, a CAD approach is advocated in three dimensions, describing high frequency effects such as current redistribution due to the skin-effect or eddy currents and the occurrence of slow-wave modes. The electromagnetic environment is described by a scalar electric potential and a magnetic vector potential. These potentials are not uniquely defined, and in order to obtain a consistent discretization scheme, a gauge-transformation field is introduced. The displacement current is taken into account to describe current redistribution and a small-signal analysis solution scheme is proposed based upon existing techniques for static fields in semiconductors. In addition methods and apparatus for refining the mesh used for numerical analysis is described.

RELATED APPLICATIONS

This application is a continuation of, and incorporates by reference inits entirety, U.S. application Ser. No. 10/630,439, filed Jul. 29, 2003,which is a continuation of U.S. application Ser. No. 09/888,868, titled“METHOD AND APPARATUS FOR SIMULATING PHYSICAL FIELDS”, filed Jun. 25,2001, now U.S. Pat. No. 6,665,849. U.S. Pat. No. 6,665,849: (1) claimsthe benefit of U.S. Provisional Application No. 60/213,764, filed Jun.23, 2000, and is hereby incorporated by reference in its entirety; and(2) is a continuation-in-part of U.S. application Ser. No. 09/328,882,titled “A METHOD FOR LOCALLY REFINING A MESH”, filed Jun. 9, 1999, nowU.S. Pat. No. 6,453,275, which claims the benefit of provisionalApplication No. 60/088,679, filed Jun. 9, 1998.

BACKGROUND OF THE INVENTION Field of the Invention

The present invention relates to a method and apparatus for simulatingfields especially electromagnetic fields, particularly useful in thecontext of analysis of interconnect structures, is presented.

Many problems in engineering, physics and chemistry require solvingsystems of partial differential equations of the type:${{{\overset{->}{\nabla}{\cdot {\overset{->}{J}}^{(k)}}} + \frac{\partial\rho^{(k)}}{\partial t}} = S^{(k)}};$k being a positive whole numberIn this equation (equation 1), J represents a flux of a substance underconsideration whose density is given by ρ and S represents some externalsource/sink of the substance. To mention a few examples:

Electrical engineering:

ρ is the charge density,

J is the current density,

S is the external charge source (recombination, generation, . . . )

Structural engineering:

Computational Fluid Dynamics:

ρ^((i))=u^((i)) the components of the fluid velocity field,

J^((i))=P−μ/3∇·u^((i)) the pressure tensor,

S^((i))=μ/ρ∇²u^((i))+F^((i))/m the external force,

∂ρ/∂t→∂ρ/∂t={right arrow over (∇)}·{right arrow over (u)} the convectivederivative

The list is not exhaustive and there exist many examples where problemsare reformulated in such a form that their appearance is as in equation(1). An important example is the Laplace equation and Poisson equation,in which{right arrow over (J)}={right arrow over (∇)}ψis the derivative of a scalar field.

There have been presented a number of methods for solving a set ofpartial differential equations as given above. All numerical methodsstart from representing the continuous problem by a discrete problem ona finite set of representative nodes in the domain where one isinterested in the solution. In other words a mesh is generated in apredetermined domain. The domain can be almost anything ranging from atleast a part of or a cross-section of a car to at least a part of or across-section of a semiconductor device. For clarification purposes, thediscussion is limited here to two-dimensional domains andtwo-dimensional meshes. This mesh comprises nodes and lines connectingthese nodes. As a result, the domain is divided into two-dimensionalelements. The shape of the elements depends amongst others on thecoordinate system that is chosen. If for example a Cartesian coordinatesystem is chosen, the two-dimensional elements are e.g. rectangles ortriangles. Using such a mesh, the domain can be introduced in a computeraided design environment for optimization purposes. Concerning the mesh,one of the issues is to perform the optimization using the appropriateamount of nodes at the appropriate location. There is a minimum amountof nodes required in order to ensure that the optimization process leadsto the right solution at least within predetermined error margins. Onthe other hand, if the total amount of nodes increases, the complexityincreases and the optimization process slows down or even can fail.Because at the start of the optimization process, the (initial) meshusually thus not comprise the appropriate amount of nodes, additionalnodes have to be created or nodes have to be removed. Adding nodes iscalled mesh refinement whereas removing nodes is called mesh coarsening.Four methods are discussed. As stated above, for clarification andsimplification purposes the ‘language’ of two dimensions is used, butall statements have a translation to three or more dimensions.

The finite-difference method is the most straightforward method forputting a set of partial differential equations on a mesh. One dividesthe coordinate axes into a set of intervals and a mesh is constructed byall coordinate points and replaces the partial derivatives by finitedifferences. The method has the advantage that it is easy to program,due to the regularity of the mesh. The disadvantage is that during meshrefinement many spurious additional nodes are generated in regions whereno mesh refinement is needed.

The finite-box method, as e.g. in A. F. Franz, G. A. Franz, S.Selberherr, C. Ringhofer and P. Markowich “Finite Boxes—A Generalizationof the Finite-Difference Method Suitable for Semiconductor DeviceSimulation” IEEE Trans. on Elec. Dev. ED-30, 1070 (1983), is animprovement of the finite-difference method, in the sense that not allmesh lines need to terminate at the domain boundary. The mesh lines mayend at a side of a mesh line such that the mesh consists of a collectionof boxes, i.e. the elements. However, numerical stability requires thatat most one mesh line may terminate at the side of a box. Therefore meshrefinement still generates a number of spurious points. The issue of thenumerical stability can be traced to the five-point finite differencerule that is furthermore exploited during the refinement.

The finite-element method is a very popular method because of its highflexibility to cover domains of arbitrary shapes with triangles. Thechoice in favor of triangles is motivated by the fact that each trianglehas three nodes and with three points one can parameterize an arbitrarylinear function of two variables, i.e. over the element the solution iswritten as:ψ(x,y)=a+b·x+c·y

In three dimensions one needs four points, i.e. the triangle becomes atetrahedron. The assembling strategy is also element by element.Sometimes for CPU time saving reasons, one performs a geometricalpreprocessing such that the assembling is done link-wise, but this doesnot effect the element-by-element discretization and assembling. Thedisadvantage is that programming requires a lot of work in order toallow for submission of arbitrary complicated domains. Furthermore,adaptive meshing is possible but obtuse triangles are easily generatedand one must include algorithms to repair these deficiencies, sincenumerical stability and numerical correctness suffers from obtusetriangles. As a consequence, mesh refinement and in particular adaptivemeshing, generates in general spurious nodes.

The finite-element method is not restricted to triangles in a plane.Rectangles (and cubes in three dimensions) have become popular. However,the trial functions are always selected in such a way that a uniquevalue is obtained on the interface. This restriction makes sense forrepresenting scalar functions ψ(x,y) on a plane.

In the box-integration method, each node is associated with an area(volume) being determined by the nodes located at the closest distancefrom this node or in other words, the closest neighbouring node in eachdirection. Next, the flux divergence equation is converted into anintegral equation and using Gauss theorem, the flux integral of thesurface of each volume is set equal to the volume integral at the righthand side of the equation, i.e. equation 1 becomes${\int_{\partial\Omega_{n}}{{\overset{->}{J}}^{(i)} \cdot {\overset{->}{\mathbb{d}}s}}} = {\int_{\Omega_{n}}{\left( {S^{(i)} - \frac{\partial\rho^{(i)}}{\partial t}} \right){\mathbb{d}^{n}x}}}$

The assembling is done node-wise, i.e. for each node the surfaceintegral is decomposed into contributions to neighboring nodes and thevolume integral at the right-hand side is approximated by the volumetimes the nodal value. The spatial discretization of the equation thenbecomes${\sum\limits_{k}{J_{lk}\frac{\partial\Omega_{lk}}{h_{lk}}}} = {\left( {S_{l}^{(i)} - \frac{\partial\rho_{l}^{(i)}}{\partial t}} \right){\Delta\Omega}_{l}}$

The advantages/disadvantages of the method are similar as for the Finiteelement method because the control volumes and the finite elements areconjugate or dual meshes. Voronoi tessellation with the Delaunayalgorithm is often exploited to generate the control volumes.

However, forming and refining the mesh is not the only problem facingthe skilled person in the solution of field theory problems. Forinstance, the on-chip interconnect structure in modern ULSI integratedcircuits is a highly complex electromagnetic system. The full structuremay connect more than one million transistors that are hosted on asilicon substrate, and containing up to seven metallization layers, andincluding interconnect splittings, curves, widenings, etc. A structureresults with a pronounced three-dimensional character. As a consequence,analytic solution methods have only limited applicability and numericalor computer-aided design methods need to be used. The continuousdown-scaling of the pitch implies that parasitic effects become a majordesign concern. Furthermore, interconnect delay will soon become themain bottleneck for increasing the operation frequencies of the fullyintegrated circuit. These observations justify an in-depth analysis ofthe interconnect problem based on the basic physical laws underlying thedescription of these systems. Whereas in the past it sufficed to extractthe parasitic behavior from the low-frequency values of thecharacteristic parameters, such as the resistance (R), capacitance (C)and inductance (L), knowledge about the modifications of theseparameters due to fast variations in time of the fields, i.e. at highfrequency, becomes mandatory. A generic method that allows one to obtainthe frequency dependence of the characteristic parameters for aninterconnect (sub-) system has been a requirement for some time.

The highest frequency in which is currently of interest is 50 GHz, whichcorresponds to a shortest wavelength of the order of one centimeter.However, this is only a current limit. For most of the interconnectswith sub-micron widths the characteristic width (length) of thestructure is therefore much smaller than the wavelengths underconsideration. In this regime one normally neglects the fulldisplacement current, but this view must be refined depending upon thematerials used [H. K. Dirks, Quasi-Stationary Fields for MicroelectronicApplications, Electrical Engineering, 79, 145-155, 1996]. Interconnectlines are typically parallel to the axes of a Cartesian grid Manhattanlike geometry. Although this is no longer true for widenings andsplittings in the lines and the vertical connections, i.e. cylindricalvias, most of the structure can be regarded in a first orderapproximation as consisting of straight orthogonal lines or bricks. Theskin effect becomes important for the upper metallization levels wherethe width of the structures is larger than the skin depth for aluminumor copper, especially at the high frequency part of the spectrum. Eddycurrents play an important role in the lossy semiconductor siliconsubstrate. It is desirable to formulate the equations for theinterconnect system in a language that is familiar to theinterconnect-designer community. In particular, variables such as thePoisson field should have their usual meaning. For time-dependent fieldsit can be achieved by selecting a specific gauge fixing. In particular,in the Coulomb gauge, the Poisson equation remains unaltered. Thenatural choice for the description of interconnect systems is the onethat uses the electric scalar potential and the magnetic vectorpotential. Small signal analysis (AC analysis) has been a successfultool for extracting compact model parameters for devices [S. E. Laux,Techniques for Small-Signal Analysis of Semiconductor devices, IEEEtrans. on computer-aided design, 4, 472-481, 1985]. Recently goodresults were obtained in using small-signal analysis [S. Jenei, privatecommunication, 2000] for the extraction of compact model parameters forthe Hasagawa system [H. Hasegawa et al. IEEE Trans. on Microwave Theoryand Techniques vol. MTT-19, 869, 1971] and similar methods are currentlyexploited for the design of spiral inductors.

Numerical analysis is well known to the skilled person, e.g. “The finiteelement method”, Zienkiewicz and Taylor, Butterworth-Heinemann, 2000 or“Numerical Analysis”, Burden and Faires, Brooks/Cole, 2001. Conventionalfinite difference numerical analysis solves three-dimensional fieldtheory problems that contain the magnetic vector potential bysuperimposing three scalar fields, representing this vector potential,whereby each scalar value is located at a node of a mesh. Finitedifference methods convert partial differential equations into algebraicequations for each node based on finite differences between a node ofinterest and a number of neighbours. These methods introduce three typesof errors. Firstly, there is the error caused by solving for a discretemesh, which is only an approximation to a continuum. The smaller themesh the higher the accuracy. Secondly, the finite difference methodsrequire an iterative solution, which is terminated after a certaintime—this implies a residual error. Thirdly, the superposition of threescalar fields is only an accurate representation of vector fields whenthe mesh size is so small that moving from one node to the next in onedirection is associated with a negligible change in the field values inthe other two dimensions. In such a case small changes of dimension inone direction may be considered as if the values of the field in theother two are constant. Where there are strongly varying fields thiscriterion can only be met where the mesh spacing is very small, i.e.there are a large number of nodes. Computational intensity increasesrapidly with the number of nodes. To a certain extent the computationalintensity can be reduced by modifying the size of mesh so that a tightmesh is only used where the divergence of the field requires this.However, varying mesh sizes places limitations on the continuity of thesolution resulting in unnecessary nodes being created to providesufficiently gradual changes. Hence, conventionally a large amount ofstorage space and high-powered computers are required to achieve anaccurate result in a reasonable amount of time.

AIM OF THE INVENTION

It is an aim of the present invention to provide numerically stablemethods and apparatus implementing these methods for simulating (i.e.calculating) field problems, e.g. electromagnetic fields.

It is a further aim of the present invention to provide numericallystable methods and apparatus implementing these methods for simulating(i.e. calculating) field problems, e.g. electromagnetic fields whichrequires less storage space and preferably less computational intensity.

SUMMARY OF THE INVENTION

The present invention provides a consistent solution scheme for solvingfield problems especially electromagnetic modeling that is based uponexisting semiconductor techniques. A key ingredient in the latter onesis the numerical solutions method based on a suitable finite differencemethod such as the Newton-Raphson technique for solving non-linearsystems. This technique requires the inversion of large sparse matrices,and of course numerical stability demands that the inverse matricesexist. In particular, the finite difference matrix, e.g. aNewton-Raphson matrix should be square and non-singular. The presentinvention provides a generic method for solving field problems, e.g.simulating electromagnetic fields, and is designed for numericalstability, in particular the solution of partial differential equationsby numerical methods.

It is an aspect of the invention that it is recognized that in order toobtain a consistent discretization scheme, meaning leading to numericalfit calculations, a dummy transformation field, also denoted gaugetransformation field or auxiliary gauge field can be introduced as adummy field and can ease computation. The dummy field can be introduceddue to the non-uniqueness of the electric and magnetic potentialsdescribing the underlying physical phenomena.

It is an aspect of the invention that it is recognized that in order toobtain a consistent discretization scheme special caution is taken inthe translation of the continuous field equations onto the discretelattice, comprising of nodes and links.

With the generic method high-frequency parasitic effects and thefrequency dependence of the characteristic parameters for aninterconnect (sub-) system can be studied but the method is not limitedthereto.

The present invention provides a method for numerical analysis of asimulation of a physical system, the physical system being describableby field equations in which a parameter is identifiable as a one-formand solving for a field equation corresponding to the parameter resultsin a singular differential operation, the method comprising:

-   -   directly solving the field equations modified by addition of a        dummy field by numerical analysis, and    -   outputting at least one parameter relating to a physical        property of the system.

The method can be formalized as follows: a method for simulating fieldsin or about a device, said method comprising the steps of:

-   -   modifying the set of field equations expressed in terms of the        vector potential of said fields to a set of modified field        equations expressed in terms of the vector potential of said        inductive fields and a dummy field; and    -   directly solving the set of modified field equations in order to        obtain the vector potential and said dummy field.

The output of the method is a field related parameter of the device,e.g. an electromagnetic parameter of the device such as a fieldstrength, a resistivity, an inductance, a magnetic field strength, anelectric field strength, an energy value. The field equations of theabove method may be the Maxwell equations. The dummy field is preferablya scalar field.

The present invention also provides a method for numerical analysis of asimulation of a physical system, the physical system being describableby Maxwell's field equations of which the following is a representation:$\begin{matrix}{{\nabla{\times \left( {\frac{1}{\mu}{\nabla{\times A}}} \right)}} = {J - {ɛ\frac{\partial}{\partial t}\left( {{\nabla V} + \frac{\partial A}{\partial t}} \right)}}} & (1) \\{{\nabla{\cdot A}} = 0} & (2) \\{{- {\nabla\left( {ɛ{\nabla V}} \right)}} = \rho} & (3) \\{E = {{- {\nabla V}} - \frac{\partial A}{\partial t}}} & (4) \\{B = {\nabla{\times A}}} & (5)\end{matrix}$

Where J and ρ are generic functions of the fields, i.e.J=J(E,B,t)  (6)ρ=ρ(E,B,t)  (7)

the method comprising:

-   -   directly solving the field equations modified by addition of a        dummy field by numerical analysis, the dummy field removing a        singularity in the numerical analysis, and    -   outputting at least one parameter relating to a physical        property of the system.

The physical property may be any of the following non-limiting list: anelectric current, a current density, a voltage difference, an electricfield value, a plot of an electric field, magnetic field value, a plotof a magnetic field, a resistance or a resistivity conductance or aconductivity, a susceptance or a suceptibility, an inductance, anadmittance, a capacitance, a charge, a charge density, an energy of anelectric or magnetic field, a permittivity, a heat energy, a noise levelinduced in any part of a device caused by electromagnetic fields, afrequency.

The above methods also include a step refining a mesh used in thenumerical analysis in accordance with an embodiment of the presentinvention.

The present invention may provide an apparatus for numerical analysis ofa simulation of a physical system, the physical system being describableby field equations in which a parameter is identifiable as a one-formand solving for a field equation corresponding to the parameter resultsin a singular differential operation, the apparatus comprising: meansfor solving by numerical analysis a modification of the field equations,the modification being an addition of a dummy field, and means foroutputting at least one parameter relating to a physical property of thesystem.

The present invention may also provide an apparatus for numericalanalysis of a simulation of a physical system, the physical system beingdescribable by Maxwell's field equations of which the following is arepresentation: $\begin{matrix}{{\nabla{\times \left( {\frac{1}{\mu}{\nabla{\times A}}} \right)}} = {J - {ɛ\frac{\partial}{\partial t}\left( {{\nabla V} + \frac{\partial A}{\partial t}} \right)}}} \\{{\nabla{\cdot A}} = 0} \\{{- {\nabla\left( {ɛ{\nabla V}} \right)}} = \rho} \\{E = {{- {\nabla V}} - \frac{\partial A}{\partial t}}} \\{{B = {\nabla{\times A}}}{where}{J = {J\left( {E,B,t} \right)}}{\rho = {\rho\left( {E,B,t} \right)}}}\end{matrix}$

the apparatus comprising:

-   -   means for directly solving the field equations modified by        addition of a dummy field by numerical analysis, the dummy field        being added to remove a singularity during the numerical        analysis, and    -   means for outputting at least one parameter relating to a        physical property of the system.

The present invention may include a data structure for use in numericalanalysis of a simulation of a physical system, the physical system beingdescribable by field equations in which a parameter is identifiable as aone-form and solving for a field equation corresponding to the parameterresults in a singular differential operation, the field equations beingmodified by addition of a dummy field, wherein the data structurecomprises the simulation of the physical system as a representation ofan n-dimensional mesh in a predetermined domain of the physical system,the mesh comprising nodes and links connecting these nodes therebydividing said domain in n-dimensional first elements whereby eachelement is defined by 2^(n) nodes, the data structure being stored in amemory and comprising representations of the nodes and the links betweennodes, the data structure also including definitions of a parameter ofthe dummy field associated with the nodes of the mesh.

A data structure for use in numerical analysis of a simulation of aphysical system, the physical system being describable by Maxwell'sfield equations of which the following is a representation:$\begin{matrix}{{\nabla{\times \left( {\frac{1}{\mu}{\nabla{\times A}}} \right)}} = {J - {ɛ\frac{\partial}{\partial t}\left( {{\nabla V} + \frac{\partial A}{\partial t}} \right)}}} \\{{\nabla{\cdot A}} = 0} \\{{- {\nabla\left( {ɛ{\nabla V}} \right)}} = \rho} \\{E = {{- {\nabla V}} - \frac{\partial A}{\partial t}}} \\{{B = {\nabla{\times A}}}{where}{J = {J\left( {E,B,t} \right)}}{\rho = {\rho\left( {E,B,t} \right)}}}\end{matrix}$

-   -   the field equations being modified by addition of a dummy field,        wherein the data structure comprises the simulation of the        physical system as a representation of an n-dimensional mesh in        a predetermined domain of the physical system, the mesh        comprising nodes and links connecting these nodes thereby        dividing said domain in n-dimensional first elements whereby        each element is defined by 2^(n) nodes, the data structure being        stored in a memory and comprising representations of the nodes        and the links between the nodes, the data structure also        including definitions of the vector potential A associated with        the links of the mesh.

The present invention also includes a program storage device readable bya machine and encoding a program of instructions for executing any ofthe methods of the present invention.

The present invention also includes a computer program product fornumerical analysis of a simulation of a physical system, the physicalsystem being describable by field equations in which a parameter isidentifiable as a one-form and solving for a field equationcorresponding to the parameter results in a singular differentialoperation, the computer program product comprising: code for solving thefield equations modified by addition of a dummy field by numericalanalysis, and code for outputting at least one parameter relating to aphysical property of the system.

The present invention also includes a computer program product fornumerical analysis of a simulation of a physical system, the physicalsystem being describable by Maxwell's field equations of which thefollowing is a representation: $\begin{matrix}{{\nabla{\times \left( {\frac{1}{\mu}{\nabla{\times A}}} \right)}} = {J - {ɛ\frac{\partial}{\partial t}\left( {{\nabla V} + \frac{\partial A}{\partial t}} \right)}}} \\{{\nabla{\cdot A}} = 0} \\{{- {\nabla\left( {ɛ{\nabla V}} \right)}} = \rho} \\{E = {{- {\nabla V}} - \frac{\partial A}{\partial t}}} \\{{B = {\nabla{\times A}}}{where}{J = {J\left( {E,B,t} \right)}}{\rho = {\rho\left( {E,B,t} \right)}}}\end{matrix}$

the computer program product comprising:

-   -   code for solving the field equations modified by addition of a        dummy field by numerical analysis, and    -   code for outputting at least one parameter relating to a        physical property of the system.

The present invention also includes a method for numerical analysis of asimulation of a physical system, comprising: transmitting from a nearlocation a description of the physical system to a remote location wherea processing engine carries out any of the method in accordance with thepresent invention, and receiving at a near location at least onephysical parameter related to the physical system. The modified fieldequations are: $\begin{matrix}{{{\nabla{\times \left( {\frac{1}{\mu}{\nabla{\times A}}} \right)}} - {\gamma{\nabla\chi}}} = {J - {ɛ\frac{\partial}{\partial t}\left( {{\nabla V} + \frac{\partial A}{\partial t} + \frac{\partial{\nabla\chi}}{\partial t}} \right)}}} & (8) \\{{{\nabla{\cdot A}} + {\nabla^{2}\chi}} = 0} & (9) \\{{- {\nabla\left( {ɛ{\nabla V}} \right)}} = \rho} & (10) \\{E = {{- {\nabla\left( {V + \frac{\partial\chi}{\partial t}} \right)}} - \frac{\partial A}{\partial t}}} & (11) \\{B = {\nabla{\times \left( {A + {\nabla\chi}} \right)}}} & (12)\end{matrix}$

-   -   where γ is non-zero a scaling factor, which guarantees matching        of dimensions.

In the present invention the introduction of the dummy field representedby χ preferably does not modify the vector potential A found from thesolution of the modified field equations when compared with the vectorpotential found from solution of the unmodified field equations. Theaccuracy of the method may be checked by comparing known algebraicsolutions of simple fields with the solution of the method according tothe present invention.

In the method the step of directly solving the set of modified fieldequations is performed by discretizing the set of modified fieldequations onto a mesh with nodes and links between said nodes. Forexample, the mesh can be a Cartesian mesh. In particular, in the methodthe vector potential is defined on the links of the mesh. The advantageof associating a field vector with the links and not with the nodesresults from the fact that links define a direction given inherently bythe form of the mesh. Hence, a field vector field is associated with anatomic vector element of the mesh. This is a more accurate simulationthan using the superposition of scalar fields to simulate a fieldvector. the present invention makes advantageous use of vector elementsin the mesh to solve the field equations more accurately. This reducesthe number of nodes required to achieve a certain accuracy. This alsomeans that the amount of memory required is reduced as well as speedingup the calculation time.

In the method the dummy field is also defined on the nodes of the meshas it is a scalar field. That is in the finite difference method thenodes are used as the reference points for values of the dummy field. Inthe method other terms in the modified field equations are expressed interms of the vector potential and the dummy field. In the method thecurl-curl operation on a vector potential on a link is expressed infunction of the vector potential on the link and the vector potentialson neighboring links of this link. The curl-curl operation on a vectorpotential on a link is expressed in function of the vector potential onthe link and the vector potentials on links, defined by wings with saidlink.

In the method the step of directly solving can exploit a Newton-Raphsonprocedure for solving nonlinear equations. In this case it is preferredto select the dummy field in order to have square non-singular matricesin the Newton-Raphson procedure.

In the method the boundary conditions may be determined by solving aMaxwell equation in a space with 1 dimension less than the space inwhich the original field equations are solved.

In a further aspect of the invention a method, i.e. the so-calledCube-Assembling Method (CAM), is disclosed for locally refining an-dimensional mesh in a predetermined domain, wherein the mesh comprisesnodes and n−1 planes connecting these nodes thereby dividing said domainin n-dimensional first elements. This method may be advantageouslycombined with other embodiments of the invention for the solution offield theory equations. The domain can be almost anything ranging fromat least a part of a car to at least a part of a semiconductor device.For clarification purposes, the present invention will be described withreference to two-dimensional domains and two-dimensional meshes but thepresent invention is not limited thereto. The shape of the elementsdepends amongst others on the coordinate system, which is chosen. Byapplying a mesh on a domain, the domain can be introduced in a computeraided design environment for optimization purposes. Concerning the mesh,one of the issues is to perform the optimization using the appropriateamount of nodes at the appropriate location. There is a minimum amountof nodes required in order to ensure that the optimization process leadsto the right solution at least within predetermined error margins. Onthe other hand, if the total amount of nodes increases, the complexityincreases and the optimization process slows down or even can fail.Because at the start of the optimization process, the (initial) meshusually thus not comprise the appropriate amount of nodes, additionalnodes have to be created or nodes have to be removed during theoptimization process. Adding nodes is called mesh refinement whereasremoving nodes is called mesh coarsening. The method of the presentinvention succeeds in adding or removing nodes locally. The assemblingis done over the elements, being e.g. squares or cubes or hypercubesdependent of the dimension of the mesh. Like the finite-box method, theCAM method is easy to program, even in higher dimensions. However, theCAM method does not suffer from the restriction that only one line mayterminate at the side of a box.

According to this aspect of the invention, a method is disclosed forlocally refining a n-dimensional mesh in a predetermined domain, whereinthe mesh comprises nodes and n−1 planes connecting these nodes therebydividing said domain in n-dimensional first elements whereby eachelement is defined by 2^(n) nodes, said method comprising at least thesteps of:

-   -   creating a first additional node inside at least one of said        first elements by completely splitting said first element in        exactly 2^(n) n-dimensional second elements in such a manner        that said first additional node forms a corner node of each of        said second elements which results in the replacement of said        first element by said 2^(n) n-dimensional second elements; and    -   creating a second additional node inside at least one of said        second elements by completely splitting said second element in        exactly 2^(n) n-dimensional third elements in such a manner that        said second additional node forms a corner node of each of said        third elements which results in the replacement of said second        element by said 2^(n) n-dimensional third elements.

In an embodiment of the invention after the mesh is locally refined,this mesh is locally coarsened.

In another embodiment of the invention, the refinement is based on anadaptive meshing strategy.

In another aspect of the invention, a program storage device isdisclosed storing instructions that when executed by a computer performthe method for locally refining a n-dimensional mesh in a predetermineddomain, wherein the mesh comprises nodes and n−1 planes connecting thesenodes thereby dividing said domain in n-dimensional first elementswhereby each element is defined by 2^(n) nodes, said method comprisingat least the steps of:

-   -   creating a first additional node inside at least one of said        first elements by completely splitting said first element in        exactly 2^(n) n-dimensional second elements in such a manner        that said first additional node forms a corner node of each of        said second elements which results in the replacement of said        first element by said 2^(n) n-dimensional second elements; and    -   creating a second additional node inside at least one of said        second elements by completely splitting said second element in        exactly 2^(n) n-dimensional third elements in such a manner that        said second additional node forms a corner node of each of said        third elements which results in the replacement of said second        element by said 2^(n) n-dimensional third elements.

In an aspect of the invention a method is disclosed for optimizing of apredetermined property of a n-dimensional structure, said methodcomprising the steps of:

-   -   creating a n-dimensional mesh on at least a part of said        structure; said mesh containing nodes and n−1 planes connecting        these nodes thereby dividing said domain in n-dimensional first        elements whereby each element is defined by 2^(n) first element;    -   refining said n-dimensional mesh by creating a first additional        node inside at least one of said first elements by completely        splitting said first element in exactly 2^(n) n-dimensional        second elements in such a manner that said first additional node        forms a corner node of each of said second elements which        results in the replacement of said first element by said 2^(n)        n-dimensional second elements;    -   further refining said n-dimensional mesh by creating a second        additional node inside at least one of said second elements by        completely splitting said second element in exactly 2^(n)        n-dimensional third elements in such a manner that said second        additional node forms a corner node of each of said third        elements which results in the replacement of said second element        by said 2^(n) n-dimensional third elements; and    -   where said n-dimensional mesh is used to create an improved        structure.

In an embodiment of the invention said structure improvements are basedon extracting said property from structure characteristics, determinedat a subset of said nodes of said mesh.

In a further embodiment of the invention said structure characteristicsare determined by solving the partial differential equations, describingthe physical behavior of said structure, on said mesh.

The present invention will now be described with reference to thefollowing drawings.

BRIEF DESCRIPTION OF THE DRAWINGS

FIG. 1 shows a schematic representation of a computing device which maybe used with the present invention.

FIG. 2 shows placement of field variables to be solved on a Cartesiangrid in accordance with an embodiment of the present invention.

FIG. 3 shows the assembly of the curl-curl-operator using 12contributions of neighboring links in accordance with an embodiment ofthe present invention.

FIG. 4 shows the assembly of the div-grad-operator using 6 contributionsof neighboring nodes in accordance with an embodiment of the presentinvention.

FIGS. 5A and 5B shows how the boundary conditions of the B-field outsidethe simulation domain is determined in accordance with an embodiment ofthe present invention.

FIG. 6 shows the numbering applied in a 2×2×2 cube case in accordancewith an embodiment of the present invention.

FIG. 7 shows the B-field of a current on a wire.

FIG. 8 shows a magnetic field around a straight conductor, as calculatednumerically (+) in accordance with an embodiment of the presentinvention, compared with the exact (−−) algebraic solution.

FIG. 9 shows how the node pointers are arranged logically in a datastructure according to an embodiment of the present invention.

FIG. 10 shows how the link pointers are arranged logically in a datastructure according to an embodiment of the present invention.

FIG. 11 shows how the cube pointers are arranged logically in a datastructure according to an embodiment of the present invention.

FIG. 12 shows the layout of a metal plug on a highly doped semiconductorused to demonstrate the methods according to the present invention.

FIG. 13 shows doping in the semiconductor region of the metal on thehighly-doped semiconductor plug.

FIGS. 14 and 15 show magnetic field plots of the static solution seen inperspective from the top and bottom plane.

FIG. 16 shows a layout of two crossing wires used to demonstrate themethods of the present invention.

FIG. 17 a layout of a square coax structure used to demonstrate themethods of the present invention.

FIG. 18 a layout of a spiral inductor structure used to demonstrate themethods of the present invention.

FIG. 19 shows the magnetic field strength in the plane of the spiralconductor of FIG. 18 as calculated by a method of the present invention.

FIGS. 20A and 20B depict the assembling strategy, according to anembodiment of the invention. The flux in link ab is composed of twoparts: a contribution from the lower rectangle (element) and acontribution from the upper rectangle (element).

FIG. 21 depicts a mesh according to an embodiment of the invention,wherein each node is associated with an area, i.e. the black area, beingdetermined by the nodes located at the closest distance from this nodeor in other words, the closest neighbouring node in each direction Eachnode is connected to at most eight different nodes in the mesh.

FIGS. 22A, 22B and 22C depict an initial mesh and this mesh after afirst and a second local refinement according to an embodiment of theinvention.

FIG. 23 depicts a transition of a mesh based on a first orthogonalcoordinate system to a mesh based on another orthogonal coordinatesystem using the method of the present invention.

FIG. 24 depicts the node balance assembling technique according to anembodiment of the invention.

FIG. 25 depicts the structure lay-out of the diode.

FIG. 26 depicts the initial square mesh of the diode.

FIG. 27 depicts the square mesh after 1 adaption sweep.

FIG. 28 depicts the square mesh after 2 adaption sweep.

FIG. 29 depicts the square mesh after 3 adaption sweep.

FIG. 30 depicts the square mesh after 4 adaption sweep.

FIG. 31 depicts the square mesh after 5 adaption sweep.

FIG. 32 depicts the square mesh after 6 adaption sweep.

FIG. 33 depicts the current-Voltage plot.

FIG. 34 is a flowchart.

FIG. 35 is a block diagram.

FIG. 36 is a flowchart.

FIG. 37 is a flowchart.

DETAILED DESCRIPTION OF ILLUSTRATIVE EMBODIMENTS

The present invention will be described with reference to certainembodiments and drawings but the present invention is not limitedthereto but only by the claims. In particular, the present inventionwill be described with reference to the solution of electromagneticfield problems especially those associated with semiconductor devicesbut the skilled person will appreciate that the present invention hasapplication to the solution of field theory problems in general and tothe solution of partial differential equations in general.

Without being limited by theory, embodiments of the present inventionrelating to the solution of electromagnetic field equations are based onthe following observations. The Maxwell equations formulated in terms ofE and B allow for a geometrical interpretation analogous to fluiddynamics. In this picture the electric field E is a one-form, in otherwords a numerical value is assigned to each path in space. The numericalvalue corresponds to the work done by the electric field when a chargewould move along the path. The magnetic field B is a two-form i.e. anumerical value is assigned to each area element that counts the numberB-field flux lines (flows) that pass through the area element.

There is a different and more abstract geometrical interpretation ofelectrodynamics. The fields E and B can be expressed as derivatives of ascalar V and vector potential A:$E = {{- {\nabla V}} - \frac{\partial A}{\partial t}}$ B = ∇×A

Now the fields E and B can be viewed as the curvature of a space. Thisis the space of phases that may be assigned to quantum fields. Thiscurvature interpretation is lacking in the older geometrical picture ofelectrodynamics.

In order to detect the strength of the electromagnetic field it sufficesto go around an infinitesimal loop and measure the mismatch between thestarting value of the phase factor and the end value of the phasefactor. In analytic calculations this interpretation has no seriousconsequences because these calculations are based comparinginfinitesimal changes in the variables going from one position toanother. Therefore, the vector potential can be regarded as a field i.e.its dependence on the space-time variables is only local. However, in acomputer calculations are made using a grid or mesh of nodes and linksand neighboring positions (grid nodes) are always a finite distanceapart. Therefore, the round trip along a closed loop for detecting theelectromagnetic field consists of line segments that are also of finitelength. The phase factor of each line segment depends on the details ofthe path and therefore the assignment of the vector potential, inaccordance with an aspect of the present invention, is done to thesepaths. In fact, the exact connection between the phase factor, the pathC and the vector potential reads as:${\varphi\left\lbrack {{{path}\quad C},A} \right\rbrack} = {\exp\frac{\mathbb{i}}{\hslash}{\int_{C}{A \cdot {\mathbb{d}x}}}}$

Only in the limit of the mesh size going to zero are there no seriousconsequences. However, the use of very small mesh sizes increases thememory requirement as well as the time for processing all these nodes ina finite difference scheme. Determining the limit as the mesh size goesto zero is not practical in numerical analysis. The present inventionpresents new solutions of the Maxwell equations arising from the newgeometrical interpretation (and any other field equations having similarsingularity problems) using the numerical analysis of finite mesh sizes.

The new geometrical interpretation of electrodynamics requiresassignment of the vector potential A to the links which are the pathsegments of the grid in the following way:A_(ij) ≈A·Δxwhere ij refers to the link between the two neighboring nodes i and j.The vector A is represented by its projections onto the three axes x, y,z and a value is assigned to each mesh link in these directions. Hence,the vectorial nature of A is maintained by assigning it to a link whichitself is a vector. The present invention therefore makes advantageoususe of the inherent vectorial nature of a grid of nodes and links innumerical analysis.

In general, if a parameter can be defined as a one-form, i.e. a mappingfrom a line segment to a number, then this parameter should berepresented in the computer code as a variable assigned to the links ofthe grid. This observation can be widened even more: If a parameter canbe identified as a two-form it should be assigned to the area elements(plaquettes) of the grid, and if a parameter can be identified as athree-form it should be assigned to volume elements of the grid. Areference providing technical background information of the abovegeometrical representations is “The Geometry of Physics”, TheodoreFrankel, Cambridge Univ. Press, 1997.

However, there remains a serious difficulty. The Maxwell equationsformulated in terms of the potentials V and A, exhibit a singularbehavior with respect to the inversion of the full differentialoperators. The fact that the differential operator is singular indicatesan underlying symmetry. This symmetry is eliminated in accordance withan aspect of the present invention by symmetry breaking conditions. Inthe present invention a method based on a dummy or auxiliary field (χ),is described to convert the singular behavior of the differentialoperator into a regular differential operator without altering thephysical content of the system of equations.

Summarizing: if a parameter can be identified as a one-form and thecorresponding coupling between nearest neighbor nodes of a mesh used fornumerical analysis results in a singularity problem, the singularity maybe alleviated by the inclusion of an auxiliary parameter withoutaltering the physical realization (solution) of the inversion problem.

FIG. 1 is a schematic representation of a computing system which can beutilized with the methods and in a system according to the presentinvention. A computer 10 is depicted which may include a video displayterminal 14, a data input means such as a keyboard 16, and a graphicuser interface indicating means such as a mouse 18. Computer 10 may beimplemented as a general purpose computer, e.g. a UNIX workstation.

Computer 10 includes a Central Processing Unit (“CPU”) 15, such as aconventional microprocessor of which a Pentium III processor supplied byIntel Corp. USA is only an example, and a number of other unitsinterconnected via system bus 22. The computer 10 includes at least onememory. Memory may include any of a variety of data storage devicesknown to the skilled person such as random-access memory (“RAM”),read-only memory (“ROM”), non-volatile read/write memory such as a harddisc as known to the skilled person. For example, computer 10 mayfurther include random-access memory (“RAM”) 24, read-only memory(“ROM”) 26, as well as an optional display adapter 27 for connectingsystem bus 22 to an optional video display terminal 14, and an optionalinput/output (I/O) adapter 29 for connecting peripheral devices (e.g.,disk and tape drives 23) to system bus 22. Video display terminal 14 canbe the visual output of computer 10, which can be any suitable displaydevice such as a CRT-based video display well-known in the art ofcomputer hardware. However, with a portable or notebook-based computer,video display terminal 14 can be replaced with a LCD-based or a gasplasma-based flat-panel display. Computer 10 further includes userinterface adapter 19 for connecting a keyboard 16, mouse 18, optionalspeaker 36, as well as allowing optional physical value inputs fromphysical value capture devices such as sensors 40 of an external system20. The sensors 40 may be any suitable sensors for capturing physicalparameters of system 20. These sensors may include any sensor forcapturing relevant physical values required for solution of the fieldproblems, e.g. temperature, pressure, fluid velocity, electric field,magnetic field, electric current, voltage. Additional or alternativesensors 41 for capturing physical parameters of an additional oralternative physical system 21 may also connected to bus 22 via acommunication adapter 39 connecting computer 10 to a data network suchas the Internet, an Intranet a Local or Wide Area network (LAN or WAN)or a CAN. This allows transmission of physical values or arepresentation of the physical system to be simulated over atelecommunications network, e.g. entering a description of a physicalsystem at a near location and transmitting it to a remote location, e.g.via the Internet, where a processor carries out a method in accordancewith the present invention and returns a parameter relating to thephysical system to a near location.

The terms “physical value capture device” or “sensor” includes deviceswhich provide values of parameters of a physical system to be simulated.Similarly, physical value capture devices or sensors may include devicesfor transmitting details of evolving physical systems. The presentinvention also includes within its scope that the relevant physicalvalues are input directly into the computer using the keyboard 16 orfrom storage devices such as 23.

A parameter control unit 37 of system 20 and/or 21 may also be connectedvia a communications adapter 38. Parameter control unit 37 may receivean output value from computer 10 running a computer program fornumerical analysis in accordance with the present invention or a valuerepresenting or derived from such an output value and may be adapted toalter a parameter of physical system 20 and/or system 21 in response toreceipt of the output value from computer 10. For example, the dimensionof one element of a semiconductor device may be altered based on theoutput, a material may be changed, e.g. from aluminium to copper, or amaterial may be modified, e.g. a different doping level in asemiconductor layer, based on the output.

Computer 10 also includes a graphical user interface that resides withinmachine-readable media to direct the operation of computer 10. Anysuitable machine-readable media may retain the graphical user interface,such as a random access memory (RAM) 24, a read-only memory (ROM) 26, amagnetic diskette, magnetic tape, or optical disk (the last three beinglocated in disk and tape drives 23). Any suitable operating system andassociated graphical user interface (e.g., Microsoft Windows) may directCPU 15. In addition, computer 10 includes a control program 51 whichresides within computer memory storage 52. Control program 51 containsinstructions that when executed on CPU 15 carry out the operationsdescribed with respect to any of the methods of the present invention.

Those skilled in the art will appreciate that the hardware representedin FIG. 1 may vary for specific applications. For example, otherperipheral devices such as optical disk media, audio adapters, or chipprogramming devices, such as PAL or EPROM programming devices well-knownin the art of computer hardware, and the like may be utilized inaddition to or in place of the hardware already described.

In the example depicted in FIG. 1, the computer program product (i.e.control program 51) can reside in computer storage 52. However, it isimportant that while the present invention has been, and will continueto be, that those skilled in the art will appreciate that the mechanismsof the present invention are capable of being distributed as a programproduct in a variety of forms, and that the present invention appliesequally regardless of the particular type of signal bearing media usedto actually carry out the distribution. Examples of computer readablesignal bearing media include: recordable type media such as floppy disksand CD ROMs and transmission type media such as digital and analoguecommunication links.

In one embodiment, the computer 10 includes certain components that cancomprise hardware, software, or a combination thereof. For example, thecomputer 10 includes a solving component for solving equations and anoutputting for outputting data.

Maxwell Equations

The interconnect modeling directly relies upon the Maxwell equations,that describe the temporal and spatial evolution of the electromagneticfields in media.

Gauss' law∇·D=ρ  (13)

Absence of magnetic monopoles∇·B=0  (14)

Maxwell-Faraday $\begin{matrix}{{\nabla{\times E}} = {- \frac{\partial B}{\partial t}}} & (15)\end{matrix}$

Maxwell-Ampere $\begin{matrix}{{\nabla{\times H}} = {J + \frac{\partial D}{\partial t}}} & (16)\end{matrix}$where D, E, B, H, J and ρ denote the electrical induction, the electricfield, the magnetic induction, the magnetic field, the current densityand the charge density, respectively.Constitutive Laws

The following constitutive equations relate the inductances to the fieldstrengths:B=μH  (17)D=εE  (18)

The constitutive equation that relates the current J to the electricfield and the current densities, is determined by the medium underconsideration. For a conductor the current J is given by Ohm's law.J=σE  (19)where the current density satisfies the current-continuity equation:$\begin{matrix}{{{\nabla{\cdot J}} + \frac{\partial\rho}{\partial t}} = 0} & (20)\end{matrix}$

In a dielectric there are no free charges. As a simplifyingapproximation, the case will be considered of a dielectric medium whoselossy effects can be neglected. In this case, no current continuityequations need to be solved in the dielectric materials and theirdielectric constants may be assumed to be real. Although this is asevere restriction, the dielectric materials that are used in back-endprocessing of semiconductor devices are sufficiently robust againstenergy absorption, in order to preserve signal integrity at thefrequencies under consideration [A. Von Hippel, Dielectric materials andapplications, Artech House, Boston, 1995]. In the semiconductingregions, the current J consists of negatively and positively chargedcarrier currents obeying the current continuity equations.$\begin{matrix}{{{\nabla{\cdot J_{n}}} - {q\frac{\partial n}{\partial t}}} = {U\left( {n,p} \right)}} & (21) \\{{{\nabla{\cdot J_{p}}} + {q\frac{\partial p}{\partial t}}} = {- {U\left( {n,p} \right)}}} & (22)\end{matrix}$

In here, the charge and current densities areρ=q(p−n+N _(D) −N _(A))  (23)J _(n) =qμ _(n) E+kTμ _(n) ∇n  (24)J _(p) =qμ _(p) pE−kTμ _(p) ∇p  (25)J=J _(n) +J _(p)  (26)and U(n,p) is the generation/recombination of charge carriers. Thecurrent continuity equations provide the solution of the variables n andp. Note that the permittivity ε in equation 18 is real whereas, for theapplications envisaged, it may be safely assumed in the following thatthe structure is non-magnetic, i.e. μ may be assumed to be equal to μ₀).Potential Description

In order to implement the equations into software algorithms, anelectric scalar potential V and a magnetic vector potential A isintroduced in the following way. From equation 14 the magnetic inductionB may be written asB=∇×(A+∇χ)  (27)where χ is an arbitrary scalar field. The presence of the field χ isclearly mathematically redundant since ∇×(∇χ)=0. Moreover, ∇χ can beabsorbed in the vector potential A.

As will demonstrated in section on the discretization scheme, the fieldχ is a key ingredient to set up a consistent discretization scheme.Inserting equation 27 into equation 15 yields: $\begin{matrix}{{{\nabla{\times \left( {E + \frac{\partial A}{\partial t} + \frac{\partial{\nabla\quad\chi}}{\partial t}} \right)}} = 0}{whence}} & (28) \\\begin{matrix}{E = {{- {\nabla\quad V}} - \frac{\partial A}{\partial t} + \frac{\partial{\nabla\quad\chi}}{\partial t}}} \\{= {{- {\nabla\left( {V + \frac{\partial\chi}{\partial t}} \right)}} - \frac{\partial A}{\partial t}}}\end{matrix} & (29)\end{matrix}$where the last equality reflects the arbitrariness in the definition ofthe scalar potential V. Insertion of equations 27 and 29 into theremaining Maxwell equations 13 and 16 gives: $\begin{matrix}{{{- \nabla} \cdot \left( {{ɛ{\nabla\quad V}} + {ɛ\frac{\partial A}{\partial t}} + {ɛ\frac{\partial{\nabla\chi}}{\partial t}}} \right)} = \rho} & (30) \\{{{\frac{1}{\mu}{\nabla{\times \left( {\nabla{\times A}} \right)}}} - {\gamma{\nabla\chi}}} = {J - {ɛ\frac{\partial}{\partial t}\left( {{\nabla\quad V} + \frac{\partial A}{\partial t} + \frac{\partial{\nabla\chi}}{\partial t}} \right)}}} & (31)\end{matrix}$where γ is a scaling factor, which should be non-zero, e.g. 1 or −1.Since A is not uniquely determined, an appropriate gauge still has to bechosen. In order to maintain a connection to the usual language andsyntax of the static modeling of interconnects, a generalized Coulombgauge such that Poisson's equation is recovered may be chosen:∇·A+∇ ²χ=0  (32)

The basic equations can now be summarized as $\begin{matrix}{{\nabla{\cdot \left( {ɛ{\nabla\quad V}} \right)}} = {- \rho}} & (33) \\{{{\frac{1}{\mu}{\nabla{\times \left( {\nabla{\times A}} \right)}}} - {\gamma{\nabla\chi}}} = {J - {ɛ\frac{\partial}{\partial t}\left( {{\nabla\quad V} + \frac{\partial A}{\partial t} + \frac{\partial{\nabla\chi}}{\partial t}} \right)}}} & (34)\end{matrix}$

It should be stressed that so far no approximations have been made. Theregular Coulomb gauge corresponds to χ=0 and is convenient for analyticcalculations. In particular, after manipulating the term ∇×∇×A as−∇²A+∇(∇·A), the last term vanishes and equation 34 becomes−∇² A=μ ₀(J+J _(D))  (35)where J_(D) is the displacement current. Analytic solution schemesaddress equation 35 as a three-fold Poisson equation. This approach isusually sustained in numerical solution schemes, distributing the threecomponents A_(x), A_(y), A_(z), over the nodes of the discrete lattice.

As indicated above there are strong arguments to associate the field Ato links. First of all, from a gauge-theoretical point of view, thefield A is the Lie algebra element that describes the phase factor of apath in real space. A successful discretization of gauge theoriesassigns the group elements, and therefore the gauge fields to links [K.G. Wilson, Confinement of Quarks, Phys. Rev. D10, 2445, 1974]. Anotherargument in favor of this association is that the vector potential canbe identified in differential geometry with a one-form, i.e. a functionon vectors, where in accordance with the present invention the vectorsare connecting two adjacent grid nodes [T. Frankel, The Geometry ofPhysics, Cambridge University Press, 1997]. With these arguments in minda gauge field variable A_(ij)=A·ê_(ij) is associated to each link whereê_(ij) is a unit vector in the direction of the link between nodes i andj.

The time evolution can be described either in real time or in theFourier domain. In one aspect of the present invention the solution tothe field equations will be performed in the Fourier domain. In order tosmooth the transition in going from the static to the dynamicdescription, a calculation scheme is provided that generates the usualcharacteristic parameters (R,C,L,G) that now become dependent on theoperation frequency ω. In the Fourier domain the potential descriptionbecomes for the selected gauge: $\begin{matrix}{{\nabla{\cdot \left( {ɛ{\nabla\quad V}} \right)}} = {- \rho}} & (36) \\{{{\frac{1}{\mu_{0}}{\nabla{\times \left( {\nabla{\times A}} \right)}}} - {\gamma{\nabla\chi}}} = {J - {j\quad{\omega ɛ}{\nabla\quad V}} + {{ɛ\omega}^{2}A} + {{ɛ\omega}^{2}{\nabla\chi}}}} & (37) \\{{{\nabla{\cdot A}} + {\nabla^{2}\chi}} = 0} & (38)\end{matrix}$Solution Scheme

Analogous to the time-dependent analysis of devices, the interconnectsystem is treated as a multi-port device with a number of ‘stand-by’operation conditions at the terminals. In particular, these conditionscan be imposed as constant voltage biases or as constant currentinjections. The stand-by conditions are assumed to be static andtherefore, firstly, the static or lowest-order solution is found. Thefrequency dependent solution is then obtained by superposition of theinput signals and the stand-by conditions. Since the magnetic field partplays an essential role in the high-frequency analysis, the magneticfield is preferably included from the start such that the appropriatedistribution of electric and magnetic energy is present in the lowestorder solution. Starting with the equations 36-38, let ω→0. This staticsolution (V₀,A₀) will correspond to the stand-by conditions. Startingwith the static solution, the different independent variables ξ(=A, V,χ, ρ, n, p) may be rewritten as a static part (with subscript index ₀)and a non-static part, denoted with a superscript hat ˆ i.e.,ξ=ξ₀+{circumflex over (ξ)}e^(jωt). Performing a Taylor series expansionand keeping only the linear terms, the result is a linearized systemthat can be solved to give the next order solution for the charge andcurrent distributions.

Static Approach

The electrostatic field, V₀, is obtained by solving the Poisson equation∇·(ε∇V ₀)=−ρ(V ₀)  (39)and the corresponding charge distribution ρ(V₀) must be calculatedself-consistently for (a) bounded surface charges on the boundarysurfaces of the dielectric regions taking into account the appropriateboundary conditions, (b) free surface charges on the boundaries of aconductor and (c) space charge in the doped semiconductor volume. Thecurrent density J₀, gives rise to the vector potential A₀, being thesolution of $\begin{matrix}{{{\frac{1}{\mu}{\nabla{\times \left( {\nabla{\times A_{0}}} \right)}}} - {\gamma{\nabla\chi_{0}}}} = {J_{0}\left( V_{0} \right)}} & (40)\end{matrix}$and submitted to the gauge condition∇·A ₀+∇²χ₀=0  (41)

For conducting media the latter equation is supplemented by∇·J ₀=0  (42)J ₀ −σE ₀=0  (43)E ₀ +∇V=0  (44)whereas in the semiconducting regions the following equations apply:ρ₀ =q(p ₀ −n ₀ +N _(D) −N _(A))  (45)J _(n0) =qμ _(n) n ₀ E ₀ +kTμ _(n) ∇n ₀  (46)J _(p0) =qμ _(p) p ₀ E ₀ −kTμ _(p) ∇p ₀  (47)∇J _(n0) =U(n ₀ ,p ₀)  (48)∇J _(p0) =−U(n ₀ ,p ₀)  (49)Linearization

In order to extract the RCLG parameters of some interconnectsub-structure, its response to a small harmonic perturbation around agiven bias operating point is considered. The bias operation point is asolution of the static set of equations. The equations that determinethe amplitudes and phases of the harmonic perturbations are obtained bylinear perturbation of the full system. Returning to equations 36-38:∇·(ε∇{circumflex over (V)})−{circumflex over (ρ)}=0  (50)1/μ∇×(∇×Â)−γ∇{circumflex over (χ)}−Ĵ+jωε∇{circumflex over (V)}−εω ² Â−εω²∇{circumflex over (χ)}=0  (51)∇·Â+∇ ²{circumflex over (χ)}=0  (52)where the sources Ĵ and {circumflex over (ρ)} must be determined by theconstitutive equations. For metals, the following equations areappropriate:∇·Ĵ+jω{circumflex over (ρ)}=0  (53)Ĵ−σÊ=0  (54)Ê+∇{circumflex over (V)}+jωÂ+jω{circumflex over (χ)}=0  (55)

For semiconductors: $\begin{matrix}{{\hat{\rho} - {q\hat{p}} + {q\hat{n}}} = 0} & (56) \\{{\hat{E} + {\nabla\hat{V}} + {{j\omega}\hat{A}} + {{j\omega}{\nabla\hat{\chi}}}} = 0} & (57) \\{{{\hat{J}}_{n} - {q\quad\mu_{n}E_{0}\hat{n}} - {q\quad\mu_{n}n_{0}\hat{E}} + {{kT}\quad\mu_{n}{\nabla\hat{n}}}} = 0} & (58) \\{{{\hat{J}}_{p} - {q\quad\mu_{p}E_{0}\hat{p}} - {q\quad\mu_{p}p_{0}\hat{E}} - {{kT}\quad\mu_{p}{\nabla\hat{p}}}} = 0} & (59) \\{{{{\nabla{\cdot {\hat{J}}_{n}}} - {{jq}\quad\omega\hat{n}} - \frac{\partial U}{\partial n}}❘_{0}{{\hat{n} - \frac{\partial U}{\partial p}}❘_{0}\hat{p}}} = 0} & (60) \\{{{{\nabla{\cdot {\hat{J}}_{p}}} - {{jq}\quad\omega\hat{p}} + \frac{\partial U}{\partial n}}❘_{0}{{\hat{n} + \frac{\partial U}{\partial p}}❘_{0}\hat{p}}} = 0} & (61)\end{matrix}$

Here, the electric field dependence of U has been suppressed but thiscan easily be taken into account. The equations describe the deviationof the system from the static solution, using the potential fields V andA, the gauge transformation field χ and the densities n en p, asindependent variables, as is illustrated in FIG. 2.

Discretization Scheme

Because of the specific geometry of the problem, the set of equations isdiscretized on a regular Cartesian grid having N nodes in eachdirection. The total number of nodes in D dimensions is M_(nodes)=N^(D).To each node may be associated D links along the positive directions,and therefore the grid has roughly D N^(D) links. This is ‘roughly’because nodes at side walls will have less contributions. In fact, thereare 2 D sides with each a number of N^((D-1)) nodes. Half the fractionof side nodes will not contribute a link in the positive direction.Therefore, the precise number of links in the lattice is M_(links)=DN^(D) (1−1/N).

As far as the description of the electromagnetic field is concerned, thecounting of unknowns for the full lattice results in M_(links) variables(A_(ij)) for the links, and M_(nodes) variables (V_(i)) for the nodes.Since each link (or node) gives rise to one equation, the naive countingis consistent. However, the gauge condition has not yet implemented. Theregular Coulomb gauge ∇·A=0, constrains the link degrees of freedom andtherefore not all link fields are independent. There are 3N³(1−1/N) linkvariables and 3N³(1−1/N)+N³ equations, including the constraints. As aconsequence, at first sight it seems that one is confronted with anoverdetermined system of equations, since each node provides an extraequation for A. However, the translation of the Maxwell-Ampere equationon the lattice leads to a singular matrix, i.e. not all rows areindependent. The rank of the corresponding matrix is 3N³(1−1/N), whereasthere are 3N³(1−1/N)+N³ rows and 3N³(1−1/N) columns. Such a situation ishighly inconvenient for solving non-linear systems of equations. Thisarises because the source terms are themselves dependent on the fields.The application of the Newton-Raphson method requires that the matricesin the Newton equation be non-singular and square. In accordance with anaspect of the present invention, the non-singular and square form of theNewton matrix can be recovered by introducing the more general gauge∇·A+∇²χ=0, where an additional field χ, i.e. one unknown per node, isincluded. Then the number of unknowns and the number of equations matchagain. In the continuum limit (N→∞), the field χ and one component of Acan be eliminated. However, on a discrete finite lattice the auxiliaryfield is essential for numerical stability. It may be concluded that thespecific gauge only serves as a tool to obtain a consistentdiscretization scheme.

It should be emphasized that the inclusion of the gauge-fixing field χshould not lead to unphysical currents. As a consequence, the χ-fieldshould be a solution of ∇χ=0.

To summarize: instead of solving the static problem∇×(∇×A)=μ₀ J(A)∇A=0  (62)the following system of equations is solved:∇×(∇×A)−γ∇χ=μ₀ J(A)∇A+∇ ²χ=0  (63)

The implementation of the gauge condition results in a unique solutionand simultaneously arrives at a system containing the same number ofequations and variables. Hence a square Newton-Raphson matrix isguaranteed while solving the full set of non-linear equations.

Differential Operators in Cartesian Grids

The div-operator integrated over a test volume ΔV_(i) surrounding a nodei can be discretized as a combination of 6 neighboring links.$\begin{matrix}{{\int_{\Delta\quad V_{i}}{{\nabla{\cdot A}}{\mathbb{d}v}}} = {\int_{\partial{({\Delta\quad V_{i}})}}{A \cdot {\mathbb{d}{\left. S \right.\sim{\overset{6}{\sum\limits_{k}}{S_{ik}A_{ik}}}}}}}} & (64)\end{matrix}$The symbol ˜ represents the conversion to the grid formulation.

The grad-operator for a link ij can be discretized as a combination of 2neighboring nodes. Integrating over a surface S_(ij) perpendicular tothe link ij gives $\begin{matrix}{\int_{\Delta\quad S_{ij}}{{{\nabla V} \cdot {\mathbb{d}{\left. S \right.\sim\frac{V_{j} - V_{i}}{h_{ij}}}}}S_{ij}}} & (65)\end{matrix}$

The grad operator for a link ij integrated along the link ij is givenby: $\begin{matrix}{{\int_{\Delta\quad L_{ij}}{{\nabla\quad V} \cdot {\mathbb{d}{\left. l \right.\sim V_{j}}}}} - V_{i}} & (66)\end{matrix}$

The curl-curl operator can be discretized for a link ij as a combinationof 12 neighboring links and the link ij itself. As indicated in FIG. 3,the field in the dual mesh, can be constructed by taking the lineintegral of the vector potential for the four ‘wings’. Integration overa surface S_(ij) perpendicular to the link ij gives $\begin{matrix}\begin{matrix}{{\frac{1}{\mu_{0}}{\int_{\Delta\quad S_{ij}}{\nabla{\times {\nabla{\times {A \cdot {\mathbb{d}S}}}}}}}} = {\frac{1}{\mu_{0}}{\int_{\partial{({\Delta\quad S_{ij}})}}{\nabla{\times {A \cdot {\mathbb{d}l}}}}}}} \\{= {\frac{1}{\mu_{0}}{\int_{\partial{({\Delta\quad S_{ij}})}}{B \cdot {\mathbb{d}l}}}}} \\{{{\sim\frac{\Lambda_{ij}}{\mu_{0}}}A_{ij}} + {\overset{12}{\sum\limits_{kl}}{\frac{\Lambda_{ij}^{kl}}{\mu_{0}}A_{kl}}}}\end{matrix} & (67)\end{matrix}$

The div-grad-operator can be discretized (FIG. 4) integrated over a testvolume ΔV_(i) surrounding a node i as a combination of 6 neighboringnodes and the node i itself. $\begin{matrix}{{\int_{\Delta\quad V_{i}}{{\nabla{\cdot \left( {ɛ{\nabla V}} \right)}}{\mathbb{d}v}}} = {\int_{\partial{({\Delta\quad V_{i}})}}{ɛ{{\nabla V} \cdot {\mathbb{d}{\left. S \right.\sim{\overset{6}{\sum\limits_{k}}{S_{ik}ɛ_{ik}\frac{V_{k} - V_{i}}{h_{ik}}}}}}}}}} & (68)\end{matrix}$Discretized Equations

The fields (V, A, χ) need to be solved throughout the simulation domain,i.e. for a semiconductor device: conductors, semiconducting regions,dielectric regions. The discretization of these equations by means ofthe box/surface-integration method gives $\begin{matrix}{{\int_{\Delta\quad S}{\left( {{\nabla{\times {\nabla{\times A}}}} - {\gamma{\nabla\chi}} - {\mu_{0}J} + {j\quad\mu_{0}{ɛ\omega}{\nabla V}} - {\mu_{0}{{ɛ\omega}^{2}\left\lbrack {A + {\nabla\chi}} \right\rbrack}}} \right) \cdot {\mathbb{d}S}}} = 0} & (69) \\{\quad{{\int_{\Delta\quad V}{\left( {{\nabla{\cdot \left( {ɛ{\nabla V}} \right)}} - \rho} \right) \cdot {\mathbb{d}v}}} = 0}} & (70) \\{\quad{{\int_{\Delta\quad V}{\left( {{\nabla{\cdot J}} + {j\omega\rho}} \right) \cdot {\mathbb{d}v}}} = 0}} & (71) \\{\quad{{\int_{\Delta\quad V}{\left( {{\nabla{\cdot A}} + {\nabla^{2}\chi}} \right) \cdot {\mathbb{d}v}}} = 0}} & (72)\end{matrix}$leading for the independent variables A, V, χ to $\begin{matrix}{{{\left( {\Lambda_{ij} - {\mu_{0}ɛ_{ij}\omega^{2}}} \right)A_{ij}} + {\overset{12}{\sum\limits_{kl}}{\Lambda_{ij}^{kl}A_{kl}}} - {\mu_{0}S_{ij}J_{ij}} + {{j\mu}_{0}{\omega ɛ}_{ij}S_{ij}\frac{V_{j} - V_{i}}{h_{ij}}} - {\left( {\gamma + {\mu_{0}ɛ_{ij}\varpi^{2}}} \right)S_{ij}\frac{\chi_{j} - \chi_{i}}{h_{ij}}}} = 0} & (73) \\{\quad{{{\overset{6}{\sum\limits_{k}}{S_{ik}ɛ_{ik}\frac{V_{k} - V_{i}}{h_{ik}}}} - Q_{i}} = 0}} & (74) \\{\quad{{{\overset{6}{\sum\limits_{k}}{S_{ik}J_{ik}}} + {{j\omega}\quad Q_{i}}} = 0}} & (75) \\{\quad{{\overset{6}{\sum\limits_{k}}{S_{ik}\left( {A_{ik} + \frac{\chi_{k} - \chi_{i}}{h_{ik}}} \right)}} = 0}} & (76)\end{matrix}$Depending on the region under consideration, the source terms(Q_(i),J_(ij)) differ.

In a conductor Ohm's law, J=σE applies, or integrated along a link ij:$\begin{matrix}{J_{ij} = {- {\sigma_{ij}\left( {\frac{V_{j} - V_{i}}{h_{ij}} + {j\omega\omega}_{ij} + {{j\omega}\quad\frac{\chi_{j} - \chi_{i}}{h_{ij}}}} \right)}}} & (77)\end{matrix}$and Q_(i) is determined by charge conservation.

For the semiconductor environment the Scharfetter-Gummel scheme can befollowed [D. L. Scharfetter, H. K. Gummel, Large scale analysis of asilicon Read diode oscillator, IEEE Trans. Elec. Devices, ED, 16, 64-77,1969]. In this approach, the diffusion equations:J=qμcE±kTμ∇c  (78)with a plus (minus) sign for negative (positive) particles areconsidered. It is assumed that the current J and vector potential A areconstant along a link and that the potential V and the gaugetransformation χ are linearly varying along the link. A local coordinateaxis u, is considered with u=0 corresponding with node i, and u=h_(ij)corresponding to node j. Integrating the equation along the link ijgives: $\begin{matrix}{J_{ij} = {{q\quad\mu_{ij}{c\left( {\frac{V_{i} - V_{j}}{h_{ij}} + {{j\omega}\left( \frac{\chi_{i} - \chi_{j}}{h_{ij}} \right)} - {{j\omega}\quad A_{ij}}} \right)}} \pm {k\quad T\quad\mu_{ij}\frac{\mathbb{d}c}{\mathbb{d}u}}}} & (79)\end{matrix}$a first-order differential equation in c, that is solved using theaforementioned boundary conditions, provides a non-linear carrierprofile. The current J_(ij) can be rewritten as $\begin{matrix}{\frac{J_{ij}}{\mu_{ij}} = {{{- \frac{\alpha}{h_{ij}}}{B\left( \frac{- \beta_{ij}}{\alpha} \right)}c_{i}} + {\frac{\alpha}{h_{ij}}{B\left( \frac{\beta_{ij}}{\alpha} \right)}c_{j}}}} & (80)\end{matrix}$using the Bernoulli function $\begin{matrix}{{{B(x)} = \frac{x}{{\mathbb{e}}^{x} - 1}}{and}} & (81) \\{\alpha = {\pm {kT}}} & (82) \\{\beta_{ij} = {q\left\lfloor {V_{i} - V_{j} + {{j\omega}\left( {ϰ_{i} - ϰ_{j}} \right)} - {{j\omega}\quad A_{ij}h_{ij}}} \right\rfloor}} & (83)\end{matrix}$The full set of equations that need to be solved is $\begin{matrix}{{- {\int_{\Delta\quad V}^{\quad}{\left\lbrack {{\nabla{\cdot J_{n}}} - {q\quad{j\omega}\quad n} - {U\left( {n,p} \right)}} \right\rbrack\quad{\mathbb{d}v}}}} = 0} & (84) \\{{\int_{\Delta\quad V}^{\quad}{\left\lbrack {{\nabla{\cdot J_{p}}} + {q\quad{j\omega}\quad p} + {U\left( {n,p} \right)}} \right\rbrack\quad{\mathbb{d}v}}} = 0} & (85)\end{matrix}$that become after discretization $\begin{matrix}{{{- {\sum\limits_{j}^{6}\quad{S_{ij}J_{nij}}}} + {q\quad{j\omega}\quad n_{i}\Delta\quad V_{i}} - {{U\left( {n_{i},p_{i}} \right)}\Delta\quad V_{i}}} = 0} & (86) \\{{{\sum\limits_{j}^{6}\quad{S_{ij}J_{pij}}} + {q\quad{j\omega}\quad p_{i}\Delta\quad V_{i}} + {{U\left( {n_{i},p_{i}} \right)}\Delta\quad V_{i}}} = 0} & (87)\end{matrix}$where all J_(ij)'s are explicitly given as a function of A_(ij), V, χ, nand p.Boundary Conditions

The simulation domain consists of an interconnect (sub-) system possiblyextended with a region of air surrounding it. Therefore, a distinctionhas to be made between boundary conditions for the simulation domain andboundary conditions for the device. For the latter, it is clear that theelectric potential V is defined on the metal terminals provided thatvoltage boundary conditions are used. The boundary conditions forsimulation domain are more subtle.

The vector potential A, needs a specific approach. It can easily be seenthat just solving equation 40 is not possible. Indeed, the left-handside of the equation is divergence-less, whereas the right hand side hasa non-vanishing divergence on the terminals, where current is enteringor leaving the structure. In order to solve this paradox, the analogoussituation of a continuity equation is considered. In the latter case,the paradox is lifted by explicitly including the external current intothe balance equation. For the curl-curl equation it is necessary toexplicitly keep track of the external B-field, i.e. by assigning toevery link at the surface of the simulation domain, a variable B_(out).At edges of the domain this field replaces two missing ‘wings’ of thecurl-curl operator, whereas at the other links of the domain surface theB-field stands for one missing ‘wing’ (FIGS. 5 a, b). The magnetic fieldoutside the simulation region B_(out) must be consistent with theexternal current distribution J_(out) over the surface of the simulationdomain. Moreover, if it is assumed that B_(out) is fully generated bythe currents that are present in the simulation problem and no externalmagnets are nearby a unique solution to equation 40 should be obtained.For this purpose, an external B_(out) perpendicular to that link isassociated with the link. Applying Ampere's law for contours asindicated in FIG. 5 a, b, the line-integral along the contour equals thetotal current crossing it, i.e. for each node $\begin{matrix}{{\frac{1}{\mu_{0}}{\oint\limits_{\partial{({\Delta\quad A_{i}})}}{B_{out} \cdot {\mathbb{d}l}}}} = {{\int\limits_{\partial{({\Delta\quad A_{i}})}}{J_{out} \cdot \quad{\mathbb{d}S}}} = I_{i}}} & (88)\end{matrix}$

Furthermore, the magnetic field must be constructed in such a way thatits divergence vanishes. For each plaquette on the simulation boundaryit implies that∇·B _(out)=0  (89)

As indicated in the case study, assembling the different matrices gives:∇·B _(out) =I ₀  (90)∇×B _(out)=0  (91)where the differential operators are acting on the link variablesB_(out) and act on the two-dimensional boundary of the simulationdomain. It should be noted that a Maxwell problem in two dimensions canbe converted into a Laplace/Poisson problem, since the vector potentialhas only one component. As a consequence, in order to solve the externalfield problem use can be made of the methods that were developed fortransmission lines. Note that the number of outside links of a regulargrid with N points in each direction, is given by M_(links) ^(out)=12(N−1)², whereas the number of nodes is given by M_(nodes) ^(out)=6 N²−12N+8 and the number of surfaces is given by M_(forces) ^(out)=6 (N−1)².This leaves M_(nodes) ^(out)+M_(forces) ^(out) equations and two more(M_(links) ^(out)) unknowns. However, the outside surface is closed andhence expressing the solenoidal character of B_(out) implies oneredundant equation. On the other hand, expressing Ampère's law forclosed paths, will result in another redundant equation, and hence aunique solution for B_(out) is obtained.

For χ it is clear from comparing equation 62 with 63 that ∇χ mustvanish. This leaves one extra degree of freedom, so that the value of χcan be chosen as equal to 0 in one point. With this choice, the valuesof χ for the other points are considered as dynamical variables, butwill result in χ=0 everywhere.

Case Studies

In order to be able to construct the differential operator matrices, achoice must be made for the numbering of nodes, edges, surfaces andvolumes. A straightforward node numbering is chosen. The numberingstarts at the corner of the box with the lowest x, y and z indices,following its neighbor nodes along the x-axis, then jumping back to thelowest x index, incrementing the y-value, and finally when the firstplane is numbered, z is incremented.

For edges, surfaces and volumes, the number associated with each numberis given by the number of the node with the smallest node index.Furthermore, S_(ij) is set to 1, and h_(ij) is set to 1, in thefollowing examples.

2×2×2 cube

A simple 8 node cube is shown in FIG. 6, where node 1 is at potentialV=1, and node 8 is at potential V=0. The following matrices representingthe differential operators can be written as $\nabla{->\begin{pmatrix}1 & {- 1} & 0 & 0 & 0 & 0 & 0 & 0 \\0 & 0 & 1 & {- 1} & 0 & 0 & 0 & 0 \\0 & 0 & 0 & 0 & 1 & {- 1} & 0 & 0 \\0 & 0 & 0 & 0 & 0 & 0 & 1 & {- 1} \\1 & 0 & {- 1} & 0 & 0 & 0 & 0 & 0 \\0 & 1 & 0 & {- 1} & 0 & 0 & 0 & 0 \\0 & 0 & 0 & 0 & 1 & 0 & {- 1} & 0 \\0 & 0 & 0 & 0 & 0 & 1 & 0 & {- 1} \\1 & 0 & 0 & 0 & {- 1} & 0 & 0 & 0 \\0 & 1 & 0 & 0 & 0 & {- 1} & 0 & 0 \\0 & 0 & 1 & 0 & 0 & 0 & {- 1} & 0 \\0 & 0 & 0 & 1 & 0 & 0 & 0 & {- 1}\end{pmatrix}}$ ${\nabla{\times {\nabla \times}}}->\begin{pmatrix}2 & {- 1} & {- 1} & 0 & {- 1} & 1 & 0 & 0 & {- 1} & 1 & 0 & 0 \\{- 1} & 2 & 0 & {- 1} & 1 & {- 1} & 0 & 0 & 0 & 0 & {- 1} & 1 \\{- 1} & 0 & 2 & {- 1} & 0 & 0 & {- 1} & 1 & 1 & {- 1} & 0 & 0 \\0 & {- 1} & {- 1} & 2 & 0 & 0 & 1 & {- 1} & 0 & 0 & 1 & {- 1} \\{- 1} & 1 & 0 & 0 & 2 & {- 1} & {- 1} & 0 & {- 1} & 0 & 1 & 0 \\1 & {- 1} & 0 & 0 & {- 1} & 2 & 0 & {- 1} & 0 & {- 1} & 0 & 1 \\0 & 0 & {- 1} & 1 & {- 1} & 0 & 2 & {- 1} & 1 & 0 & {- 1} & 0 \\0 & 0 & 1 & {- 1} & 0 & {- 1} & {- 1} & 2 & 0 & 1 & 0 & {- 1} \\{- 1} & 0 & 1 & 0 & {- 1} & 0 & 1 & 0 & 2 & {- 1} & {- 1} & 0 \\1 & 0 & {- 1} & 0 & 0 & {- 1} & 0 & 1 & {- 1} & 2 & 0 & {- 1} \\0 & {- 1} & 0 & 1 & 1 & 0 & {- 1} & 0 & {- 1} & 0 & 2 & {- 1} \\0 & 1 & 0 & {- 1} & 0 & 1 & 0 & {- 1} & 0 & {- 1} & {- 1} & 2\end{pmatrix}$ ${\nabla{\cdot \nabla}}->\begin{pmatrix}3 & {- 1} & {- 1} & 0 & {- 1} & 0 & 0 & 0 \\{- 1} & 3 & 0 & {- 1} & 0 & {- 1} & 0 & 0 \\{- 1} & 0 & 3 & {- 1} & 0 & 0 & {- 1} & 0 \\0 & {- 1} & {- 1} & 3 & 0 & 0 & 0 & {- 1} \\{- 1} & 0 & 0 & 0 & 3 & {- 1} & {- 1} & 0 \\0 & {- 1} & 0 & 0 & {- 1} & 3 & 0 & {- 1} \\0 & 0 & {- 1} & 0 & {- 1} & 0 & 3 & {- 1} \\0 & 0 & 0 & {- 1} & 0 & {- 1} & {- 1} & 3\end{pmatrix}$ $\nabla{->\begin{pmatrix}{- 1} & 0 & 0 & 0 & {- 1} & 0 & 0 & 0 & {- 1} & 0 & 0 & 0 \\1 & 0 & 0 & 0 & 0 & {- 1} & 0 & 0 & 0 & {- 1} & 0 & 0 \\0 & {- 1} & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & {- 1} & 0 \\0 & 1 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & {- 1} \\0 & 0 & {- 1} & 0 & 0 & 0 & {- 1} & 0 & 1 & 0 & 0 & 0 \\0 & 0 & 1 & 0 & 0 & 0 & 0 & {- 1} & 0 & 1 & 0 & 0 \\0 & 0 & 0 & {- 1} & 0 & 0 & 1 & 0 & 0 & 0 & 1 & 0 \\0 & 0 & 0 & 1 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 1\end{pmatrix}}$Note that the diagonal elements of the curl-curl-operator, have thevalue 2, because from the 4 ‘wings’ of FIG. 3, only 2 remain (the twoother are outside the simulation domain). The current can be found bysolving equations 42-44. In order to find expressions for B_(out),Ampère's equation is used for the contour integrals for B_(out). Fornode 1 for instance the result is that B_(out) ⁽¹⁾+B_(out) ⁽⁵⁾+B_(out)⁽⁹⁾=I⁽¹⁾. The same kind of equations for all nodes gives formally ∇.B_(out)=I. Gauss' law for the outside magnetic field becomes for thefront surface of the cube in FIG. 6, B_(out) ⁽¹⁾−B_(out) ⁽²⁾−B_(out)⁽⁵⁾+B_(out) ⁽⁶⁾=0 or formally for all surfaces: ∇×B_(out)=0. This systemof equations results in a unique solution for B_(out).

16×16×2 cube

In order to check the physical consequences of the way the method dealswith boundary conditions, a 16×16×2-node cube is simulated in which aconductor box of one volume element is implemented. At four nodes at oneside of the conductor the voltage V=1 is applied, while at the otherside the voltage is kept V=0. Hence a current will flow, characterizedby the solution of equations 42-44 in the conducting area. This solutionof J₀ determines J_(out) and B_(out) can be found as a solution ofequations 90-91. Next it is possible to calculate the magnetic vectorpotential in the simulation domain by solving equation 40 (FIG. 7). Themagnetic field in the simulation domain which is expected to change as1/r, (r representing the distance to the conductor center). A goodcorrespondence with the theory is recovered as shown in FIG. 8.

The skilled person will appreciate certain aspects of the aboveembodiments of the present invention. A method is provided for thedescription and the analysis of the electromagnetic behavior of on-chipinterconnect structures, by using small-signal analysis. This avoidssolving Helmholtz equations and still gives us information on thestructure to describe effects like the current redistributions, theimpact of high frequencies on the characteristic parameters, slow wavemodes etc. A formulation of the Maxwell equations is used that is basedon a potential approach. Furthermore, the potential fields are assignedto links. This approach has severe consequences for solving the fieldequations, which are resolved by the method in accordance with thepresent invention. In order to solve the Maxwell equations a square andnon-singular Newton-Raphson matrix is needed, and such matrices areprovided by including an extra dummy potential field χ. This fieldhowever will not change the physical solution. In fact, a dedicatedgauge fixing procedure has been presented to accommodate for thenumerical stability. The magnetic vector potential is calculated bysolving a curl-curl operator equation. This task has been carried out ina box-like example of a current carrying wire. The simulated magneticfield shows the behavior of a magnetic field generated by a wire anddemonstrates the consistency and correctness of the proposed method.

Another interesting result concerns the boundary conditions. Theinclusion of the latter is dictated by the conversion of the continuumequations to the discretized equations. Consistency of the boundaryconditions demands that a separate Maxwell problem be solved in a domainof dimension D−1=2.

One aspect of the present invention is the efficient use of memoryspace. In accordance with an embodiment of the present invention datastructures are created in a memory of a computer system which areclosely associated with the numerical analysis methods described above.One possible implementation of such data structures is given below. Thedata structures are a representation of a mesh having links connectingnodes in a mesh structure. The implementation is based on a mesh formedby cubes but the present invention is not limited to cubes. Theimplementation makes use of pointers however the present invention isnot limited to pointer based systems but may include any method ofreferring to other memory locations.

node (see FIG. 9)

This structure is used to stock nodes in a list. struct node { double x; double y : double z ; struct node *next ; unsigned int nPointers ;unsigned int number ; } ;

All properties of the nodes (x, y, z, . . . ) are stored in thisstructure. The properties can be: V, the Poisson potential, ρ the chargedensity, N the dopant concentration, n and p the electron and holeconcentration, T the temperature and χ the dummy field. In accordancewith an aspect of the present invention the zero-forms and three-formsare associated with the nodes of the mesh. The nodes are internallyplaced in a linked list, where each node points to the next node and thelast node points to NULL.

Fields x Position on the X-axis. y Position on the Y-axis. z Position onthe Z-axis. next Pointer to the next node in the list. nPointers Numberof cubes that point to this node. number The nodenumber.link (FIG. 10)

This structure is used to stock links in a list. structlink { structnode *node1; struct node *node2;// Pointer to node2 structlink *link1;// Pointer to child link1 struct link *link2;// Pointer tochild link2 struct link *next;// Pointer to next link unsigned intnPointers;// Number of cubes that point to this link. Unsigned intnumber; };

All properties of the links are stored in this structure. In particularvalues for those elements of the fields such as the vector potential Awhich are associated with links are stored in this structure. Theproperties can be: A the magnetic vector potential, J, the currentdensity (carrier density in semiconductor substrates), E the electricfield. The links are identified by 2 nodes. In accordance with an aspectof the present invention, the one-forms are associated with the links.The links are internally placed in a linked list, where each link pointsto the next link and the last link points to NULL. A link can have 2childLinks. The pointers link1 and link2 point to these children. Ifthese pointers are NULL, the link doesn't have any children.

Fields node1 Pointer to the first node of a link. node2 Pointer to thesecond node of a link. link1 Pointer to the first of the childLinks.link2 Pointer to the second of the childLinks. next Pointer to the nextnode in the list. nPointers Number of cubes that point to this link.number The linknumber.cube (see FIG. 11)

This structure is used to stock cubes in a list. structcube { unsignedint number; struct cube *cube[8]; struct node *node[8]; struct link*link[12]; struct cube *next; struct cube *parent; };

The links are identified by number. This is the number of the cube.Internally, the cubes are organized in several linked lists. All cubeswith the same generation are stored in the same linked list.

Fields number The cubenumber. cube [8] An array of 8 pointers tochildCubes. Either a cube has eight childs, or it has none. If thepointers are set to NULL, the cube doesn't have childs. node [8] Anarray of 8 pointers to the nodes of a cube. Every cube has 8 nodes, toidentify it's boundaries. link [12] An array of pointers to the 12 mainlinks of a cube. Since links can have 2 children, a cube can have morethan 12 indirect links. next Pointer to the next cube in the list. If itis NULL, this is the last cube in the list for this generation. parentPointer to the parent of the cube. If the pointer is set to NULL, thisis the “biggest cube”, and doesn't have a parent, since it is no child.cubeListPointer

This structure is used to point to a cubeList. structcubeListPointer {cube *cube; struct cubeListPointer *next;

Internally, cubes are organised in several linked lists. All cubes withthe same generation are stored in the same linked list. Something isrequired to point to all these lists. This is what a cubeListPointerdoes, it points to a cubeList and to the next cubeListPointer. So thefirst cubeListPointer points to the cubeList generation 1, the next tothe list with generation 2, . . . and so on.

Fields cube A pointer to the first cube in a cubeList. If it is set toNULL, there is no list appended to the cubeListPointer yet. next Apointer to the next cubeListPointer. If it's NULL, this is the lastcubeListPointer in the list.lastNumbers

This structure is used to keep track of the last nodenumber, linknumberand cubenumber. struct lastNumbers { unsigned int lastNodeNr; unsignedint lastLinkNr; unsigned int lastCubeNr ; } ;

The last nodenumber, linknumber and cubenumber are stored here, so thatnodes/links/cubes can easily be appended and thenodenumber/linknumber/cubenumber can be filled in.

Fields lastNodeNr The highest nodenumber at a certain moment. lastLinkNrThe highest linknumber at a certain moment. lastCubeNr The highestcubenumber at a certain moment.

The calculation method for the pointers is given below.

In the following a detailed description of practical applications of themethods of the present invention are described.

In order to extract the RCLG parameters of an interconnectsub-structure, its response to a small harmonic perturbation around agiven bias operating point is considered, which is a solution of thestatic set of equations. The equations that determine the amplitudes andphases of the harmonic perturbations are obtained as linearperturbations of the full system. Returning to equations (36-38), oneobtains∇(ε∇V _(R) −εωA ₁−εω∇χ₁)+ρ_(R)=0  (92)∇(ε∇V ₁ −εωA _(R)−εω∇χ_(R))+ρ₁=0  (93)∇×∇×A _(R)−μ₀εω² A _(R)−μ₀ J _(R)−μ₀ εω∇V ₁−(γ+μ₀εω²)∇χ_(R)=0  (94)∇×∇×A ₁−μ₀εω² A ₁−μ₀ J ₁+μ₀ εω∇V _(R)−(γ+μ₀εω²)∇χ₁=0  (95)∇²χ_(R) +∇A _(R)=0  (96)∇²χ₁+∇A₁=0  (97)where the sources J_(R), J_(I), ρ_(R) and ρ₁ must be determined by thenon-linear constitutive equations.Boundary Conditions

Continuing the discussion on boundary conditions above, the vectorpotential A, needs a specific approach. It can easily be seen that justsolving equation 40 is not possible. Indeed, the left hand side of theequation is ∇²χ, whereas the right hand side has a non-vanishingdivergence on the terminals, where current is entering or leaving thestructure. However, as was argued above the dummy field χ=0 is part ofthe solution. In order to solve this paradox the analogous situation ofa continuity equation is considered. In the latter case, the paradox islifted by explicitly including the external current into the balanceequation.

The external currents, impinging perpendicular on the boundary surfaceof the simulation domain, carry their own circular magnetic field. Sucha magnetic field is described by a component of the vector potentialparallel to the impinging current. Therefore the boundary condition forthe vector potential will be equal to zero for all links that are in theboundary surface, ∂Ω, of the simulation domain, Ω, whereas linkspointing orthogonally inwards from the enclosing surface are part of theset of the unknown variables that should be solve,A_(ij)=0 i,jε∂Ω  (98)

The boundary condition for the χ-field will be that χ=0 at the enclosingsurface. The Laplace problem on a closed surface with these boundaryconditions guarantees that χ=0, everywhere.

The boundary conditions for the scalar potential V are a mixture ofDirichlet and Neumann boundary conditions. At the contacts Dirichletboundary conditions are assumed whereas at the remaining part of theenclosing surface Neumann boundary conditions are assumed. Thisassumption implies that no perpendicular electric field exists for theseparts of ∂Ω. If a contact is placed on a semi-conducting region, it isassumed that this contact is also ohmic. Therefore, the boundarycondition for a semiconductor contact isφ_(n)|_(c)=φ_(p)|_(c) =V| _(c)  (99)where φ_(p) and φ_(n) are the quasi-fermi level for the hole andelectron concentrations at the contact. Furthermore, it is assumed thatcharge neutrality holds at the contact, i.e. p−n+N=0 and np=n_(i) ².Strictly speaking, these assumptions are valid only for contactsattached to highly-doped regions, otherwise one would have to deal withSchottky contacts. However, within the framework of back-end structuresimulations, this assumption is valid since the contacts tosemiconducting regions usually are connected to highly-doped source ordrain regions or polysilicon gates.Interface Conditions

In general, the structures consist of insulating, semiconducting andmetallic regions. As a consequence, there will be four types ofinterface nodes, i.e.

insulator/metal interface nodes

insulator/semiconductor interface nodes

metal/semiconductor interface nodes

insulator/semiconductor/metal ‘triple’ points

At the metal/semiconductor interface nodes, the idealized interfaceSchottky contact condition is implemented, as for a boundary conditionfor a semiconductor region, by setting φ_(p)=φ_(n)=V_(metal), whereV_(metal) is the value of the Poisson potential at the metal side of theinterface. The Poisson potential at the semiconductor side of theinterface is V_(semi)=V_(metal)−δV, where δV represents the contactpotential between the two materials. Using $\begin{matrix}{{n = {n_{i}\exp\frac{q}{kT}\left( {V - \phi_{p}} \right)}}{p = {n_{i}\exp\frac{q}{kT}\left( {\phi_{p} - V} \right)}}} & (100)\end{matrix}$and applying the neutrality condition p−n−N=0, where N=N_(D)−N_(A) forp-type semiconductor regions (N<0) results in: $\begin{matrix}{{\delta\quad V} = {\log\quad\left( {{- \frac{N}{2n_{i}}}\left( {1 + \sqrt{1 + \frac{4n_{i}^{2}}{N^{2}}}} \right)} \right)}} & (101)\end{matrix}$and for n-type semiconductor regions (N>0) $\begin{matrix}{{\delta\quad V} = {{- \log}\quad\left( {\frac{N}{2n_{i}}\left( {1 + \sqrt{1 + \frac{4n_{i}^{2}}{N^{2}}}} \right)} \right)}} & (102)\end{matrix}$

At a metal/semiconductor interface node there is one variable(V_(metal)) that needs to be solved. The equation for this variableassigned to the node i is the current-continuity equation,$\begin{matrix}{{\sum\limits_{j}^{\quad}\quad{J_{ij}S_{ij}}} = 0} & (103)\end{matrix}$where J_(ij) is the current density in discretized form for the link(ij) from node i to node j, and S_(ij) the perpendicular cross sectionof the link (ij). Note that for an idealized Schottky contact thePoisson potential is double-valued.

At metal/insulator interface nodes continuity of the Poisson potentialis assumed. For these nodes there is, apart from the variables A and χ,one unknown V_(i), and the corresponding equation is thecurrent-continuity equation. The Poisson equation determines theinterface charge, ρ_(i) and can be obtained by post-processing, once Vis determined.

At insulator/semiconductor interface nodes there are three unknowns tobe determined, V, n and p. These variables are treated in the usual wayas is done in device simulation tools, i.e. the Poisson equation issolved self-consistently with the current-continuity equations for n andp, while V is continuous at the insulator/semiconductor interface.

At triple point nodes, the Poisson potential is triple-valued. Onearrives at different values depending on the material in which oneapproaches the node. For computational convenience the value of thePoisson potential in the insulator at the midpoint between V_(metal) andV_(semi) is taken, i.e. $\begin{matrix}{{\lim\limits_{x->x_{tr}}\quad{V_{insul}(x)}} = {{V_{metal}\left( x_{tr} \right)} - {\frac{1}{2}\delta\quad V}}} & (104)\end{matrix}$

The interface conditions for the vector field A and the ghost field χare straightforward. The choice of the gauge condition, equation 38, isindependent of specific material parameters. During assembling ofequation 37, the current associated to each link is uniquely determinedin an earlier iteration of the Gummel loop and therefore for the vectorpotential (as well as χ) is single-valued.

Scaling Considerations

In order to program the equations on a computer, an appropriate scalingmust be performed. A generalization of the ‘de Mari’ scaling, asdescribed in A. de Mari, An Accurate Numerical Steady-State OneDimensional Solution of the PN-Junction, Solid-State Electronics, 11,33-58, 1968 is adopted. Let λ be the scaling parameter for lengths,n_(i) the scaling parameter for doping and carrier concentrations andlet the thermal voltage V_(T)=[kT/q] act as a scaling parameter for thePoisson field and the fermi levels. Then from Poisson's equation oneobtains that $\begin{matrix}{\lambda = \sqrt{\frac{ɛ_{0}V_{T}}{{qn}_{i}}}} & (105)\end{matrix}$

From the constitutive equations for the carrier currents the scalingfactor for the mobility, σ_(μ) is obtained: $\begin{matrix}{\sigma_{\mu} = \frac{\sigma_{D}}{V_{T}}} & (106)\end{matrix}$

There is still the freedom to set one scaling parameter. The scalingparameter for the diffusion constant, σ_(D)=1 [(m²)/sec] is fixed. Thenthe scaling parameter for the time, τ, is given by $\begin{matrix}{\tau = \frac{\lambda^{2}}{\sigma_{D}}} & (107)\end{matrix}$

Furthermore, from the scaling factor for the diffusivity the scalingfactor for the current density is obtained: $\begin{matrix}{\sigma_{J} = \frac{{qn}_{i}\sigma_{D}}{\lambda}} & (108)\end{matrix}$

The frequency scaling factor is inverse to the time scaling factor, i.e.σ_(ω)=[1/(τ)]. The scaling factor for A and χ follows from thegeneralized formula for the electric field, $\begin{matrix}{{\sigma_{A} = \frac{\tau\quad V_{T}}{\lambda}}{\sigma_{\chi} = {\tau\quad V_{T}}}} & (109)\end{matrix}$

The scaling of the curl-curl equation for A leads to the dimensionlessconstant, K=[1/(c²)] ([(λ)/(τ)])² where c is the speed of light invacuum. From table I, K is a rather small number that makes it suitablefor using it as a perturbation expansion parameter. TABLE I Generalized‘de Mari’ scaling factors. Variable Name Value Unit Temperature T 300 KPoisson field V_(T) 2.5852 × 10⁻²  V Concentration n_(i) 10¹⁶ m⁻³ Lengthλ 1.1952 × 10⁻⁵  m Diffusion σ_(D) 1 [(m²)/sec] constant Mobility σ_(μ)38.6815 [(m²)/V sec] Current density σ_(J) 134.0431 [C/(m² sec)] Time τ1.4286 × 10⁻¹⁰ sec Electric field σ_(E) 2162.8670 [V/m] Frequency σ_(ω)6.9994 × 10⁹  sec⁻¹ Conductance σ_(σ) 6.1974 × 10⁻²  [C/Vmsec] Velocityσ_(v) 8.3662 × 10⁴  [m/sec] Chi-scaling σ_(χ) 3.6934 × 10⁻¹² V secA-scaling σ_(A) 3.0900 × 10⁻⁷  [V sec/m] K-factor K 7.7879 × 10⁻⁸ dimensionless

After scaling the curl-curl equation takes the following form∇×∇×{right arrow over (A)}−γ∇χ=K({right arrow over (J)}−jε _(r) ω∇V+ε_(r)ω² A+ε _(r)ω²∇χ)  (110)The constant γ can be used as a tuning parameter to improve convergenceof the linear solvers. It should be noted that γ may not be zero, sincefor γ=0, the equations for A and χ decouple such that the matrix blockfor A becomes singular.Numbering Schemes

The novelty of the new approach for solving the equations for the scalarand the vector potentials is the association of the vector potentialvariables to the links or connections of the discretization grid. Thisrequires that not only the grid nodes receive a unique pointer, but alsothe grid links. In order to become familiar with this new situation, amethod for assigning unique pointers to the grid nodes and the gridlinks in Cartesian grids is presented.

For a Cartesian grid with N_(node)=k_(x)×k_(y)×k_(z) nodes, a uniquenode pointer is obtained byL _(node) =n _(x)+(n _(y)−1)×k _(x)+(n _(z)−1)×k _(x) ×k _(y)  (111)where n=(n_(x),n_(y),n_(z))εIN³ points to a specific node in the gridand 1≦n_(x)≦k_(x), 1≦n_(y)≦k_(y) and 1≦n_(z)≦k_(z).

Given a node pointer L_(node), the vector n can be reconstructed usingthe following algorithm if (L_(node) = N_(node)) then n_(x) =k_(x), n_(y) = k_(y), n_(z) = k_(z) else J₀ = INT((L_(node) − 1)/(k_(x)× k_(y))) n_(z) = J₀−1 K₀ = L_(node)−J₀ ×k_(x) ×k_(y) L₀ =INT((K₀−1)/k_(x)) n_(y) = L₀+1 n_(x) = K₀ − L₀ × k_(x) endif

In a similar way a unique pointer, L_(link) can be assigned to eachlink. Given the node n and a direction i=1, 2, 3, a link pointer can beobtained by the following algorithm, since each link is based in somenode n, and points in a given positive direction, i. if (i = 1) then  if(n_(x) = k_(x)) then illegal input  else   L_(link) = n_(x) + (k_(x)−1)× (n_(y)−1) + (k_(x)−1) × k_(y) × (n_(z)−1)  endif elseif (i = 2) then if (n_(y) = k_(y)) then illegal input  else   L_(link) = (k_(x)−1) ×k_(y) × k_(z) + n_(x) + (k_(x)) ×(n_(y)−1) + k_(x) ×(k_(y)−1) ×(n_(z)−1) endif elseif (i = 3) then  if (n_(z) = k_(z)) then illegal input  else  L_(link) = 2×k_(x)×k_(y)×kz − (k_(y) +k_(x))×k_(z) + n_(x) +k_(x)×(n_(y)−1)+ k_(x)×k_(y)×(n_(z)−1)  endif endifUsing the INT-function, (n,i) can be extracted from L_(link), as wasdone for the nodes.Solver Requirements

Since a large number of equations need to be solved simultaneously, i.e.Poisson's equation, the current-continuity equations, the curl-curlequation and the equation for χ, both the static, real and imaginaryparts, a Gummel iterative procedure is followed for solving this system.In particular, since the frequency, ω is fixed, the problem is stillthree-dimensional. Whereas the Poisson's equation and thecurrent-continuity equations can be treated similarly as in devicesimulation programs such as D. L. Scharfetter, H. K. Gummel, Large ScaleAnalysis of a Silicon Read Diode Oscillator, IEEE Trans. Elec. Devices,ED, 16, 64-77, 1969 and E. M. Buturla, P. E. Cottrell, B. M. Grossmanand K. A. Salsburg, Finite-Element Analysis of Semiconductor Devices:the FIELDAY program, IBM Journal on Research and Development, 25,218-231, 1981, the equations for A and X require a different handling.Furthermore, the largest system that needs to be solved is also the pairof equations for (A, χ). These equations can not be loaded iterativelyinto the Gummel flow because the ∇x∇x operator would lead to a singularalgebraic system. However, the simultaneous assembling with the equationfor χ, results in an algebraic system having a regular matrix that canbe inverted.

After testing a number of different linear solvers with and withoutpre-conditioning, it turned out that the most robust method was thesymmetric successive-over-relaxation (SSOR) pre-conditioner, takingω_(SOR)=1.2, combined with the conjugate-gradient squared (CGS)iterative solver. The parameter γ in equation 110 was taken equal toone.

The memory requirements for the sparse storage of the Newton-Raphsonmatrix can be obtained as follows. Suppose there areN_(node)=k_(x)×k_(y)×k_(z) nodes andN_(link)=3×k_(x)×k_(y)×k_(z)−(k_(x)×k_(y)+k_(x)×k_(z)+k_(y)×k_(z))links. Each (interior) link interacts with 12 neighboring links and 2nodes in the assembling of the curl-curl equation. This implies that aeach link generates 1+12+2=15 non-zero entries in the Newton-Raphsonmatrix. The scalar equation for the χ-variable in each node interactswith 6 χ's in the neighboring nodes and with 6 link variables sited atthe links connecting the node to the nearest-neighbor points. Therefore,each χ field induces 1+6+6=13 non-zero entries in the Newton-Raphsonmatrix. The total number of non-zero entries can be estimated (ignoringsurface subtractions) as N_(non) _(—) _(zero)=15×N_(link)+13×N_(node).Table II gives a few numerical examples. TABLE II Storage requirementsfor the Newton-Raphson matrix. k_(x) k_(y) k_(z) N_(nodes) N_(links)N_(non)_zero 10 10 10 1000 2700 49060 10 10 100 10.000 27900 516340 10100 100 100.000 288.000 5.430.520 100 100 100 1000.000 2.970.00057.073.600

EXAMPLES

A number of examples are presented demonstrating that the proposedpotential formulation in terms of the Poisson field V, the vector fieldA and the dummy field χ, is a viable method to solve the Maxwell fieldproblem. All subtleties related to that formulation, i.e. thepositioning of the vector potential on links, and the introduction ofthe ghost field χ, have already been described above in constructing thesolutions of the static equations. Therefore, a series of examples inthe static limit are presented.

Metal Plug on Highly-Doped Silicon

The first example concerns the electromagnetic behavior of a metal plugon (highly-doped) silicon. This example addresses all subtleties thatare related to metal-semiconductor, metal-insulator andsemiconductor-insulator interfaces as well as triple lines. Thesimulation region (10×10×10 μm³) consists of two layers. A layer of thesilicon (5 μm) is highly doped at the top, using a square mask of 4×4μm² at the center. Above the silicon there is 5 μm oxide with a metalplug of 4×4 μm².

In FIG. 12, the structure is sketched. A Gaussian doping profile isbelow the metal plug and the concentration (at the surface of thesimulation domain) is plotted in FIG. 13. The voltage drop over the plugis 0.2 Volts. The resistivity of the metal is taken 10⁻⁸ Ωm. In FIGS. 14and 15, the magnetic field is presented. Whereas the metal plug carriesthe current in the top layer, it is observed that the field has astrength decaying as ˜1/r. In the bottom layer the current spreads outand this leads to a flat B-field intensity. In the table III the resultsare listed of some characteristic parameters of the simulation.

The energies have been calculated in two different ways and goodagreement is observed. This confirms that the new methods underlying thefield solver are trustworthy. The χ-field is zero within the numericalaccuracy, i.e. χ˜O (10⁻¹⁴). TABLE III Some characteristic results forthe metal/semiconductor plug. ELECTRIC ENERGY$\frac{1}{2}ɛ_{0}{\int_{\Omega}\quad{\mathbb{d}{vE}^{2}}}$ 2.41890E−17 J$\frac{1}{2}{\int_{\Omega}\quad{{\mathbb{d}{v\rho}}\phi}}$ 2.55777E−17 JMAGNETIC ENERGY$\frac{1}{2\mu_{0}}{\int_{\Omega}\quad{\mathbb{d}{vB}^{2}}}$ 4.84510E−23J $\frac{1}{2}{\int_{\Omega}\quad{{\mathbb{d}{vJ}}.A}}$ 4.90004E−23 JRESISTANCE Resistance 1.49831E+4 ΩCrossing Wires

The second example concerns two crossing wires. This example addressesthe three-dimensional aspects of the solver. The structure is depictedin FIG. 16 and has four ports. In the simulation one port is placed at0.1 volt and the other ports are kept at zero volt. The current is 4Ampere. The simulation domain is 10×10×14 μm³. The metal lines have aperpendicular cross section of 2×2 μm². The resistivity is 10⁻⁸ Ωm. Intable IV, some typical results are presented. TABLE IV Somecharacteristic results for two crossing wires. ELECTRIC ENERGY$\frac{1}{2}ɛ_{0}{\int_{\Omega}\quad{\mathbb{d}{vE}^{2}}}$ 1.03984E−18 J$\frac{1}{2}{\int_{\Omega}\quad{{\mathbb{d}{v\rho}}\phi}}$ 1.08573E−18 JMAGNETIC ENERGY$\frac{1}{2\mu_{0}}{\int_{\Omega}\quad{\mathbb{d}{vB}^{2}}}$ 2.89503E−11J $\frac{1}{2}{\int_{\Omega}\quad{{\mathbb{d}{vJ}}.A}}$ 2.82824E−11 JRESISTANCE Resistance 0.25 ΩSquare Coaxial Cable

To show that also inductance calculations are adequately addressed, theinductance per unit length (L) is calculated of a square coaxial cableas depicted in FIG. 17. The inductance of such a system with innerdimension a and outer dimension b, was calculated by: $\begin{matrix}{{l \times \frac{1}{2}{LI}^{2}} = {{\frac{1}{2\mu_{0}}{\int_{\Omega}^{\quad}{B^{2}\quad{\mathbb{d}v}}}} = {\frac{1}{2}{\int_{\Omega}^{\quad}\quad{{\mathbb{d}{vJ}} \cdot A}}}}} & (112)\end{matrix}$

with l denoting the length of the cable. As expected, for large valuesof the ratio r=b/a, the numerical result for the square cable approachesthe analytical result for a cylindrical cable, L=[(μ₀ ln(b/a))/(2π)].TABLE V Some characteristic results for a square coaxial cable. L L a b(cylindrical) (square) μm μm b/a (nH) (nH) 2 6 3 220 255 1 5 5 322 329 17 7 389 390 1 10 10 461 458Spiral Inductor

A spiral inductor, as shown in FIG. 18 was simulated. This structurealso addresses the three dimensional aspects of the solver. Thecross-section of the different lines is 1 μm×1 μm. The overall size ofthe structure is 8 μm×8 μm and the simulation domain is 23×20×9 μm³.

In FIG. 19, the intensity of the magnetic field is shown at height 4.5μm. TABLE VI Some characteristic results for the spiral inductor.ELECTRIC ENERGY$\frac{1}{2}ɛ_{0}{\int_{\Omega}\quad{\mathbb{d}{vE}^{2}}}$ 2.2202E−18 J$\frac{1}{2}{\int_{\Omega}\quad{{\mathbb{d}{v\rho}}\phi}}$ 2.3538E−18 JMAGNETIC ENERGY$\frac{1}{2\mu_{0}}{\int_{\Omega}\quad{\mathbb{d}{vB}^{2}}}$ 3.8077E−13J $\frac{1}{2}{\int_{\Omega}\quad{{\mathbb{d}{vJ}}.A}}$ 3.9072E−13 JRESISTANCE Resistance 0.54 Ω

The above three-dimensional field solution method has been programmed insoftware for on-chip passive structures. A TCAD software environment wasbuilt such that arbitrary Manhattan-like structures can be loaded,calculated and the results can be visualized. The multi-layer stack of aback-end process results in different material interfaces. The inclusionof the interface conditions in the TCAD software was done such that theelectromagnetic response is accurately simulated. The treatment of thedomain boundaries was done based on the idea that energy should notenter or leave the simulation domain except for the contact ports. Thenumerical implementation requires that all variables are scaled and thede Mari scaling was extended to include the vector field A and the ghostfield χ. Numbering schemes for the nodes and the links were given andthe solver requirements have been specified.

In calculating the above examples it is wasteful of nodes to have thesame mesh size over the whole of the area/volume of interest. In certainareas/volumes, e.g. where the field intensities change rapidly or inareas/volumes of great importance it is preferred to have a smaller meshsize. However, in other areas it is economical in memory size andcomputing time to have a wider mesh spacing.

In a further embodiment of the invention a method is disclosed forrefining a mesh. This method may be combined advantageously with theprevious field calculation methods or may be used in the calculation offield problems by other methods, e.g. in fluid dynamics, mechanics etc.This embodiment is therefore not limited in its use to the above filedcalculation methods. As an example the application of the method to arectangular 2-dimensional mesh in a predetermined domain will bedescribed but the present invention is not limited to this number ofdimensions. The rectangular mesh comprises nodes and one-dimensionalplanes, i.e. lines called links, connecting these nodes. As a result,the domain is divided into 2-dimensional rectangular first elementswhereby each element is defined by 2² nodes, i.e. rectangles. Theassembling is performed over these rectangles which makes sense despitethe fact that the rectangles have four nodes. Indeed, if the solutionwere modeled as a linear function over a rectangle, the number ofconditions for finding the coefficients would result into anoverdetermined problem. Therefore, the interpolation is dependent on thepurpose for which the interpolation is used. For post-processing one maycomplete the mesh with a Delaunay tessellation and exploit linearinterpolation, however, for solving the system of equations, only theend point values of each link or side enter the equations. A link (seeFIG. 20) is a connection between two adjacent nodes and forms a side ofa rectangle. In fact such a link is shared by a first and a secondrectangle, said first rectangle and said second rectangle being adjacentrectangles, and has a two-fold purpose. The link serves as the fluxcarrier for the first as well as the second rectangle, e.g. the toprectangle and the bottom rectangle. However, the flux of said firstrectangle and the flux of said second rectangle are considered to benon-interacting and therefore, they satisfy the superposition principle.This is illustrated in FIG. 20. It is this observation that allows for astable and correct assembling strategy using rectangles instead oftriangles. In this assembling strategy, particularly in two dimensions,each node is connected to at most eight different sites in the lattice,since each rectangle can generate a connection to two different nodes.These nodes may all be different, although this is not a necessity.

The basic features of a CAM algorithm in accordance an embodiment of thepresent invention are introduced with an example. Suppose one wants tosimulate the electrical potential W(x) in a rectangular grid definingsome device, given that the electrical system is described by thePoisson equation with D(x)=εE(x) the electrical field and ρ(x) theelectrical charge source and x the place coordinate. The rectangulardevice is represented by a domain D (FIG. 24). A first relation betweenthe electrical field and the electrical potential is given. A secondrelation between the electrical charge source and the electricalpotential is given. As such the Poisson equation can be expressed interms of the electrical potential.∇·D(x)=ρ(x)

∇² W(x)=ρ(W(x))

For said simulation the above equation is discretized on a rectangulargrid or mesh. Said grid divides the domain D in a set of smaller domainsD_(i). The union of said domains D_(i) (D₁−D₁₃) gives D. The equation isnow written for each node of the grid. This is further illustrated forthe node 9, central in the domain D. Around said node, four rectanglesD₅, D₇, D₆ and D₄ are recognized. In each of said rectangles areas A₁,A₄, A₃ and A₂ (shaded in FIG. 24) are defined by taking the middle ofthe links, defined between the central node and the appropriate nodes ofthe rectangles. These nodes are denoted in FIG. 24 as black colourednodes and numbered. Eight such nodes are defined. Said areas togetherdefine an overall area A=A₁+A₂+A₃+A₄. The above equation is nowintegrated over said area A and Stokes theorem is applied. Theright-hand side is approximated by the electrical charge on said centralnode. The left-hand side is replaced by a contour-integral of theelectrical field along the boundaries of the area A. Said boundariescomprises of twelve lines of length d_(i) (d₁-d₁₂). Only the electricalfield along the links to the central node can be used. In the inventedmethod the electrical field contributions along each such link areseparated in two contributions. For instance the electrical fieldcontribution along the horizontal link at the left side of the centralnode is defined to comprise of a first contribution I₄₉, running fromnode 4 to the central node and a second contribution I₈₉, running fromnode 8 to the central node. I₄₉ is related to the rectangle D₆ while I₈₉relates to D₄. Said contributions are multiplied with the line length ofthe appropriate line, being the lines, orthogonal to the electricalfield contribution under consideration. For the central node thisresults in the following equation, denoted the node balance:I ₄₉ ×d ₁ +I ₉₃ ×d ₂ +I ₉₇ ×d ₄ +I ₉₆ ×d ₅ +I ₉₂ ×d ₇ +I ₁₉ ×d ₈ +I ₅₉×d ₁₀ +I ₈₉ ×d ₁₁ −ρ ₉ ×A=0

For each of the nodes of the mesh such an equation is written. In eachequation the relation between the electrical charge ρ_(i) and theelectrical field I_(ij) with the electrical potential is introduced. Assuch a set of nonlinear equations (ρ_(i) will depend itself nonlinearlyon W_(i)) in the variables W_(i), being the electrical potential at eachnode of the mesh, is obtained. Said set of equations is then solved byusing an iterative procedure such as a Newton-Raphson scheme. In moregeneral terms I_(ij) is denoted a link-current and ρ_(i)×A is denotedsource contribution.

The CAM algorithm can be written schematically in the following form:program flux_solver call setup rectangle_init an initial mesh ofrectangles is generated 1 call solve_on_rectangles the equations aresolved on the current mesh call refine_rectangles mesh refinement usingthe CAM method if (refinement_need) go to 1 the refinement is repeatedtill a refinement criterion is met predetermined  stop

The function solve_on_rectangles can be written schematically in thefollowing form:  function solve_on_rectangles   do for each rectangle   find_variable_in_nodes    assemble link_current    do for each node    assign link_current to node_balance     assign source_contributionto node_balance    enddo   enddo solve_equations

Note that other orderings of the do-loops are also possible. Thesolve_equations routine can be any nonlinear equation solver.

According to this embodiment of the invention, a method is disclosed forlocally refining a rectangular 2-dimensional mesh in a predetermineddomain, wherein the mesh comprises nodes and lines connecting thesenodes thereby dividing said domain in 2-dimensional first elementswhereby each element is defined by 4 nodes. Particularly, this method,can locally refine an initial mesh comprising 2-dimensional firstelements, i.e. rectangles. This method comprises at least the steps of:

-   -   creating a first additional node inside at least one of said        first rectangles by completely splitting said first rectangle in        exactly four second rectangles in such a manner that said first        additional node forms a corner node of each of said second        rectangles which results in the replacement of said first        rectangle by said four second rectangles; and    -   creating a second additional node inside at least one of said        second rectangles by completely splitting said second rectangle        in exactly four third rectangles in such a manner that said        second additional node forms a corner node of each of said third        rectangles which results in the replacement of said second        rectangle by said four third rectangle.

This first additional node is created somewhere inside a rectangle. Thiscan be anywhere inside a rectangle and thus not on a link (side) of therectangle. Particularly, this first additional node can be created inthe center of the rectangle. Regardless of the exact location of thefirst additional node, this first rectangle is split completely inexactly four new rectangles, i.e. second rectangles. In other words, thesum of the areas of these four second rectangles completely coincideswith the area of this first rectangle and therefore this first rectanglecan be deleted. This also holds in an analogous way for the secondadditional node.

In another embodiment the method of the present invention is embedded inan adaptive meshing strategy. Adaptive meshing is straightforward forn-dimensional meshes comprising n-dimensional elements whereby eachelement is defined by 2^(n) nodes. Particularly adaptive meshing can beeasily implemented for a two-dimensional mesh comprising rectangles.There are no restrictions on the number of links ending on another linkas e.g. in the finite-box method. Therefore, extra algorithms forsmoothing the mesh after the adaptive meshing are not required becausethere is no metastasis of spurious nodes into irrelevant regions; inother words the refinement remains locally. In FIG. 22, a possibleresult of such a meshing strategy is shown after two cycles ofadaptation.

An important problem in simulation is the use of different physicalmodels on different length scales, or alternatively, using the samephysical models at different scales but with modified model parameters.In the latter case, the model parameters follow renormalization groupflows, as in K. Wilson “Confinement of Quarks” Phys. Rev. D10, 2445(1974). Both scenarios can be applied in restricted domains, withoutblurred transition regions, because the refinement scheme does notgenerate extra nodes which necessitate subsequent smoothing steps.

Concerning the issue of numerical stability of the method of the presentinvention, it is demonstrated that sets of equations of the type${{{\overset{->}{\nabla}{\cdot \overset{->}{J^{(k)}}}} + \frac{\partial\rho^{(k)}}{\partial t}} = S^{(k)}};$k being a positive whole numbercan be solved on a mesh obtained after a series of iterations. The keyobservation is that each mesh resembles a Kirchhoff network.Particularly, when considering a 2-dimensional rectangular mesh, eachcurrent flows along a side (link) of an element and each sideaccumulates contributions from two adjacent elements (rectangles)sharing this particular link. Moreover, the currents arising from thesetwo elements do not interact and therefore one may assemble each elementindependently. One has to keep in mind that the expressions fordiscretized currents in terms of the end point field values generatesemi-definite Newton matrices. An adaptive meshing algorithm accordingto the method of the present invention, i.e. based on therenormalization group refinement method, is created and tested on aseries of devices. No signals revealing instability are detected. Themost comprehensive simulation is performed on a MOSFET, using thehydrodynamic model for holes and electrons. Solving five equations withthe help of a simultaneous Newton solver for a mesh comprising 2000nodes, no convergence slow-down was observed. The results agreed within1% with results obtained with conventional schemes.General Orthogonal Coordinate Systems

In three dimensions there exists 12 different classes of orthogonalcoordinate systems. The orthogonality allows for applying thesubdivision of any cell in 4 subcells by taking the mean of the minimumand maximum of the coordinate variables. The orthogonality is asufficient reason for the guaranty that one does not have to includesmoothing nodes. Another interesting application concerns the transitionfrom one coordinate system to another. When using the method of thepresent invention, there are no restrictions on the number of linksending on another link. Therefore, the absence of smoothing nodesimplies that only transition nodes have to be generated at the preselected interval region. This is done by detecting the nodes incoordinate system A and in the transition region and by assigningrefinement nodes to the coordinate system B. In the same way the nodesof coordinate system B lying in the transition region become new nodesfor coordinate system A. This is illustrated in FIG. 23. Aninterpolation function needs to be chosen in order to avoidoverdetermination of the system of nodal equations. In other words: Thecoordinate systems A and B are matched by transition boundaryconditions.

In another embodiment of the present invention, a mesh is created usingthe locally refinement method of the present invention. Thereafter themesh can be locally coarsened. When solving systems of partialdifferential equations of the type${{{\overset{->}{\nabla}{\cdot \overset{->}{J^{(k)}}}} + \frac{\partial\rho^{(k)}}{\partial t}} = S^{(k)}};$k being a positive whole numberit is often observed that nodes are generated in the adaptive meshingstrategy, during the iteration process towards the finally solution,whose existence becomes obsolete, when arriving at the final solution.These nodes are artifacts of the solution procedure and not relevant fordescribing the final solution with sufficient degree of accuracy. Withthe method presented here it becomes rather straightforward to clean upthe final mesh from these nodes. For this task it is necessary to assignto each node its hierarchy during the generation of the mesh. Theinitial mesh comprises nodes and first elements and all these nodesobtain generation index 0, i.e. so-called parent nodes. The additionalfirst nodes, which are obtained in the first generation, are given theindex 1. These are children nodes. The process is continued by thegeneration of grandchildren and after X iterations, nodes from the X-thgeneration are obtained. Mesh coarsening can now be easily executed bydeleting the nodes from the last generation and check for sustainedaccuracy. All parent nodes of a child node, which has been deleted, maynow be checked for obsoleteness, and so on.

The Cube-Assembling Meshing method is now demonstrated with thesimulation of a microelectronic structure. A detailed exposure of theadaptive meshing strategy in a realistic application is presented. Theselected structure is an implanted diode with side contact. The purposeof this example is two-fold: (1) to demonstrate that the new methodgives acceptable results, (2) to demonstrate that the method has beenimplemented in a two-dimensional device simulator, which has served as aresearch tool for predicting highly unconventional devices, such ashetero-junction type vertical transistors, multi layer HEMTS as well asmore conventional structures, such as CMOS devices in operation pointswhere effects are present which are very difficult to simulate, such ascarrier heating.

All these structures have in common that commercial software tools donot provide sufficient reliable data for the parameters which are ofinterest to the process engineer. This situation is due to the fact thatsoftware development progresses in parallel with technology development.Having available a source code for the simulation of the internaldynamics inside devices, an ongoing activity has been to improve thealgorithms which are exploited in solving the equations underlying thedevice dynamics. Although many code improvements deal with a gradualextension towards more accurate models to be implemented, occasionallydramatic jumps in code improvements are realized by implementing piecesof code which involve a new understanding of the mathematical machinerywhich applied.

In the past decade it is observed three instances of major breakthroughevents in solving the semi-conductor device equations on the computer.In 1988, a method for simulating abrupt hetero structures, byprogramming finite jumps in the matching conditions of the nodal valuesof the scalar functions in the finite-element method. In 1993, a methodwas constructed for obtaining smooth solutions as well as a robustiterative scheme for finding the solutions of the energy-balanceequations by realizing that CGS solvers require a semi-definiteNewton-Raphson Jacobian matrix. The third ‘quantum leap’ in improvingthe algorithm for finding the solutions of the semiconductor deviceequations was invented in 1998.

With decreasing device dimensions, there is an increasing need for theinclusion of quantum effects in the simulation method. The laws ofquantum-mechanics are of a substantially different character as the lawsof classical physics, e.g. even in first order quantization, wave-likeequations need to be incorporated, which are very different from Poissonequation and balance equations in general which are of the diffusivetype. This situation leads to reconsideration of the discretizationschemes, with having in mind a method which decouples the macroscopicphysics from the microscopic physics. The underlying idea is thatquantum effects are important in particular regions of the device butthat other parts the device could be described by the classical methods.Furthermore the quantum regime should be entered by locally refining themesh in order to take into account the fact that the quantum effect aredominate on a small distance scale. A renormalization grouptransformation should realize the connection between the macroscopic andmicroscopic domains. Unfortunately, local mesh refinements in adaptivemeshing schemes are always accompanied by spurious nodes which take careof a smooth transition between the fine and the coarse mesh. Suchspurious nodes are lethal to the developed ideas for incorporating thequantum physics in the simulation method. Therefore, the question,whether it would be possible to incorporate a discretization scheme forthe classical device equations, which may be submitted to localrefinement without generating a metastasis of the spurious refinementnodes, has raised and been solved by the development of the CAM-scheme.The method combines the finite-element method with the box-integrationmethod and is applicable in an arbitrary number of dimensions. Themethod is limited to meshes consisting of the union of patches oforthogonal coordinate frames.

The device under consideration consists of a np-diode, whose thirddimension is much larger as the first and the second dimension.Therefore, a two-dimensional cross section suffices for thedetermination of the diode characteristics. In a perpendicular plane,the diode has dimensions 10 times 5 μm squared. A p-type doped region ofsize 5 times 2.5 μm squared is allocated in the upper left corner with acontact named top, and the remaining region is n-typed doped. There is aside-wall contact named side connected to the n-type region. When thediode is forward biased, a current flow is expected from the side-wallcontact to the top contact. This current is not expected to beone-dimensional due to the perpendicular orientation of the contacts.The adaptive meshing algorithm is expected to allocate refinement nodessuch that the current paths can be accurately traced. In FIG. 25 thedevice lay-out is presented as it is drawn by the PRISM structuregenerator. The PRISM structure generator is a pre-processor which allowsthe user to design the geometric aspects as well as other issues, suchas material selection, contact positioning and interface positioning.

The input file diode.dat for the simulator PRISM is presented below.MESH diode.str SQUARE *===========================doping-profiles============== UNIFORM−1.00E19 0.0 2.5001 5.0 5.0 UNIFORM 1.00E19 0.0 0.0 5.0 2.5 UNIFORM1.00E19 5.0001 0.0 10.0 5.0 *===================================================== DIEL_CON sili 111.9 DIEL_CON oxid 1 3.8 * intrinsic concentration of Si: INTR_CAR sili1 1.3E10 TEMP 26.0 MOB sili 11 400. 1500.0 GF sili 03 f RECOM sili 15000.0 5000.0 BIAS top 0.0 side 0.0 ANAL BIHT CPUTIM 600.0 * max nonewton loops LOOPS 100 $ solver accuracy: ACC 1.0d-14 1.0d-14 1.0d-301.0d-30 1 500 1.0d-30 1.0d-35 NORM 1.0d-3 1.0d-2 5.0d-2 1.0d-2 5.0d-2PRINT diode.out 3 PLOT doping 6 doping SOLVE

At the first run of a new structure the output of the PRISM structuregenerator is loaded into the PRISM simulation. This is the filediode.str as indicated in the the second line. The file diode.str isprinted below. In the section STRUCTURE-LINES for each line the fourthparameter gives the number of line divisions. With this parameter acrude initial mesh for performing initial calculations towards the finalsolution are generated. These number can be provided interactively withthe structure generator. <PRIMI> <TITLE> <none> <GRID>  1.000000e+015.000000e+00 2.500000e−01 2.500000e−01 <GEOMETRY> 1  <Points>  0.000000e+00 1.250000e+00 1   1.000000e+01 1.250000e+00 1  2.500000e+00 5.000000e+00 1   2.500000e+00 0.000000e+00 1  0.000000e+00 2.500000e+00 1   1.000000e+01 2.500000e+00 1  5.000000e+00 5.000000e+00 1   5.000000e+00 0.000000e+00 1  0.000000e+00 0.000000e+00 1   1.000000e+01 0.000000e+00 1  1.000000e+01 5.000000e+00 1   0.000000e+00 5.000000e+00 1  <Lines>  0.000000e+00 1.250000e+00 0.000000e+00 2.500000e+00 1   1.000000e+011.250000e+00 1.000000e+01 0.000000e+00 1   2.500000e+00 5.000000e+005.000000e+00 5.000000e+00 1   2.500000e+00 0.000000e+00 0.000000e+000.000000e+00 1   0.000000e+00 2.500000e+00 0.000000e+00 5.000000e+00 1  1.000000e+01 2.500000e+00 1.000000e+01 1.250000e+00 1   5.000000e+005.000000e+00 1.000000e+01 5.000000e+00 1   5.000000e+00 0.000000e+002.500000e+00 0.000000e+00 1   1.000000e+01 0.000000e+00 5.000000e+000.000000e+00 1   1.000000e+01 5.000000e+00 1.000000e+01 2.500000e+00 1  0.000000e+00 5.000000e+00 2.500000e+00 5.000000e+00 1   0.000000e+000.000000e+00 0.000000e+00 1.250000e+00 1 <EXTEND> 0.000000e+00<STRUCTURE>   <Points> 16   1 2.500000e+00 5.000000e+00   2 0.000000e+002.500000e+00   3 2.500000e+00 2.500000e+00   ....   ....   <Lines> 24  1 1 3 5 0 0 0   2 2 3 5 0 0 0   3 2 4 5 0 0 0   ....   ....  <Areas> 9  1 3 4 1 2 1 1   2 1 5 6 7 1 1   ....   ....  <MIC>  <Materials> 1   1sili  <Contacts> 2   7 top   9 side   <Interfaces> 0 <END>

The execution of PRISM with the input file diode.dat generates two filesrelated to the mesh construction. The first file diode.sqr is a detaileddescription which contains all structure information for applying thecube assemble method. Part of the file diode.sqr is printed below. TITLE<none> NODE 256  2.50000000E+00 5.00000000E+00  .00000000E+002.50000000E+00  2.50000000E+00 2.50000000E+00  .00000000E+005.00000000E+00  .....  ..... SQUARES 225   19  131  132  20  1  1  212 216  50  49  1  1  128  20  1  32  1  1  245  249  250  246  1  1 .....  ..... INTERFACE 0 CONTACT 17  4  1  1  1  7  2  ..  ..  ..  ..LOG_NAMES  2  0  1 C 1 1top C 2 1side M 1 1sili END

Furthermore, there is printed a file diode.gpt, which is suitable fordrawing the initial mesh using xgnuplot. The plotting is presented inFIG. 26. In order to set up the initial square mesh from the structurefile the initial run of the simulation is done without putting bias onthe contacts. The BIAS card in the diode.dat file takes care of thisaspect. A run without bias usually is fast and is suitable for a syntaxanalysis and name matching of the input files diode.dat and diode.str.Errors are reported in the output file diode.err. The initial solutionshould also be quickly obtained (i.e. a limited number of iterationsshould be used) otherwise there is likely a problem with the dopingconditions. Part of the error file diode.err is presented below. SectionTITL  has been read in. Section NODE  has been read in. Section ELEM has been read in. Section CONT  has been read in. Section LOG_(—)  hasbeen read in. CONTACT information:  contact # 1 : length = .2500E+01 um contact # 2 : length = .2500E+01 um PRISM Version 4.0 rev 0, File runof: 1-Dec-98 15:43 <none> DIEL_CON SILI 1 11.9 DIEL_CON OXID 1 3.8INTR_CAR sili 1 1.3E10 TEMP 26.0 MOB SILI 11 400. 1500.0 Energy balancemobility model E-> T (Si) GF SILI 03 f  No surface scattering included Field dependent mobility: Arora RECOM SILI 1 5000.0 5000.0 Shockley-Read-Hall Recombination ANAL BIHT CPUTIM 600.0 LOOPS 100 ACC1.0d-14 1.0d-14 1.0d-30 1.0d-30 1 500 1.0d-30 1.0d-35 NORM 1.0d-3 1.0d-25.0d-2 1.0d-2 5.0d-2 PRINT diode.out 3 DATA contact nos units  1  2contact name  top  side This mesh contains: 256 nodes, 225 squares. L2norms:- Poisson 2.5127E+03 RESOUT: #loops, ERR1, ERR3, ALPHA, BETA,<r0,C.p(k)> || 0 |2.51E+03| .00E+00||  .00E+00 .00E+00 .00E+00|| || 2|7.01E−20| .00E+00||  1.00E+00 4.55E−20 7.01E−20||   dpsi  3.70E−04 atnode 37 Time is  .02 UPDATE1: loop iij= 0 damp TK= .10000E+01 Time is .02 applied bias V  .00000E+00  .00000E+00 L2 norms:- Poisson 4.1071E−01 RESOUT: #loops, ERR1, ERR3, ALPHA, BETA, <r0,C.p(k)> || 0|4.11E−01|  .00E+00|| 1.00E+00 4.55E−20 7.01E−20|| || 2 |2.41E−27| .00E+00|| 1.00E+00 4.57E−20 2.41E−27||   dpsi  6.85E−08 at node 37 Timeis  .02 UPDATE1: loop iij= 0 damp TK= .10000E+01 Time is  .02 THE CGS-PCMETHOD STOPPED DYNAMICALLY. ALFA= −.244083E−30 NO OF ITERATIONS = 18DAMP1: damping factors TK,TK1 = .10000E+01 .10000E+01 UPDATE dpsi=1.57E−10 d0p= 1.11E−10 d0n= 1.87E−14  at nodes( 40) ( 71) ( 113) L2norms:- Poisson 2.1668E−08 L2 norms: Hole eqn 1.5673E−08 L2 norms: Holetemp 1.5482E−07 L2 norms: Elec eqn 2.4257E−08 L2 norms: Elec temp2.2549E−07   Total 2.7590E−07   Time is  .07 CONVERGED: Currentimbalance 7.0425D-12 ICLU: maximum |dii|_LLU =  .58987E+05 minimum|dii|_LU =  .96854E−05 ICLU: maximum |dii|_AN =  .10000E+01 minimum|dii|_AN =  .16953E−04   at A row  3  at A row  6 ICLU: Number ofcorrected entries = 4 RESOUT: #loops, ERR1, ERR3, ALPHA, BETA,<r0,C.p(k)> || 0 | 5.38E−10 | 6.47E−10 || −5.29E−02 −1.28E+00−8.27E−27|| THE CGS-PC METHOD STOPPED DYNAMICALLY.ALFA= −.534617E−30 NOOF ITERATIONS = 31 DAMP1: damping factors TK,TK1 = .10000E+01 .10000E+01UPDATE dpsi= 1.25E−10 d0p= 1.28E−10 d0n= 1.04E−10   at nodes( 69) ( 6) (39)    dtemp= 1.44E−17 dtemn= 1.94E−17   at nodes( 178) ( 38) 16 16Creating file doping.ddca Creating file doping.ddca New mesh: Ny = 16 Nx= 16 with  0 interpolated values

A summary file of the results is also produced which provides the biasesat the contact and the resulting currents. The file is printed below.PRISM Version 4.0 rev 0, File run of: 1-Dec-98 15:43  DIEL_CON SILI 111.9  DIEL_CON OXID 1 3.8  INTR_CAR sili 1 1.3E10  TEMP 26.0  MOB SILI11 400. 1500.0  Energy balance mobility model E-> T (Si)  GF SILI 03 f No surface scattering included  Field dependent mobility: Arora  RECOMSILI 1 5000.0 5000.0  Shockley-Read-Hall Recombination  ANAL BIHT CPUTIM 600.0  LOOPS 100  ACC 1.0d-14 1.0d-14 1.0d-30 1.0d-30 1 5001.0d-30 1.0d-35  NORM 1.0d-3 1.0d-2 5.0d-2 1.0d-2 5.0d-2  PRINTdiode.out 3  DATA  contact nos units  1  2  contact name  top  side CONVERGED: Current imbalance 9.5643D-15  CONVERGED: Current imbalance7.0425D-12  applied bias V  .00000E+00 .00000E+00  Sheet charge Ccm−1 −6.19394E−23 1.32711E−22  Hole current Acm−1  6.26567E−10 5.57167E−31 Elec current Acm−1  1.72112E−27 −6.19524E−10  Integrated electron conc 3.7624E+12  cm−1  Integrated hole conc  1.2374E+12  cm−1  Total runtime> .917E−01mins

Having obtained an initial mesh as well as a zero bias solution, one mayproceed along different tracks. First a mesh refinement according todoping criteria and/or electric field criteria may be done. Othercriteria do not make much sense at zero bias. One may also first ramp upthe bias to some particular value and apply mesh adaption at the finalstage. Which method is selected is not relevant for demonstrating thefeasibility of the Cube-Assembling Method, since in both cases themethod should final mesh should result into a L stable and convergentsolution scheme. Which new nodes are ultimately participating depends onthe history of the application of the refinement criteria and the orderof applying first adaption and then refinement or first ramping and thenadaption may effect the presence of final nodes.

Following the first option, i.e. applying adaption on the zero-biassolution, refinement according to the doping and electric-field profileis obtained. In FIG. 27 to FIG. 32 six succesive meshes are obtained byresubmitting the zero-bias solution to the refinement tests based on thedoping profile. In the following table the number of new nodes is given.Cycle New nodes Total 0 0 256 1 78 334 2 153 487 3 306 793 4 286 1079 56 1085 6 7 1092

A current-voltage plot is presented in FIG. 33. The convergence speed issimilar to the zero-bias case, which demonstrates that thecube-assembling method works properly. As such, it was demonstrated thata new assembling strategy based on a synthesis of the finite-elementmethod and the box-integration method which allows mesh refinementwithout spurious nodes, does give a algorithm which is stable, robustand convergent.

While the invention has been shown and described with reference topreferred embodiments, it will be understood by those skilled in the artthat various changes or modifications in form and detail may be madewithout departing from the scope and spirit of this invention. Forexample, although the embodiments of the present invention have beendescribed generally with reference to Cartesian grids, the presentinvention may be applied to any form of grid used in numerical analysis.Further, although the present invention has been described withreference to the numerical analysis of Maxwell's equations it appliesequally well to the solution of partial differential equations,including the refinement of the mesh used in such solutions.

1. A method of numerical analysis of a simulation of a physical system,the physical system being describable by field equations in which aparameter is identifiable as a one-form and solving for a field equationcorresponding to the parameter results in a singular differentialoperation, the method comprising: directly solving the field equationsmodified by addition of a dummy field by numerical analysis, andoutputting at least one parameter relating to a physical property of thesystem.